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SYSTEM AND METHOD FOR PERFOSMmC A 
THREE-DIMENSIONAL VIRTUAL EXAMINATION OF OBJECTS, SUCH 

AS INTERNAL ORGANS 

SPECIFICATION 

5 TECHNICAL FIELD 

The present invention relates to a system and method for performing a 
volume based three-dimensional virtual examination, and more particularly relates to 
a system which offers enhanced visualization and navigation properties. 

BACKGROUND OF THE INVENTION 
1 0 Colon cancer continues to be a maj or cause of death throughout the 

world. Early detection of cancerous growths, which in the human colon initially 
manifest themselves as polyps, can greatiy improve a patient's chance of recovery. 
Presently, there are two conventional ways of detecting polyps or otiier masses in the 
colon of a patient. The first method is a colonoscopy procedure, which uses a flexible 
1 5 fiber-optic tube called a colonoscope to visually examine the colon by way of physical 
rectal entty with the scope. The doctor can manipulate the tube to search for any 
abnormal growths in the colon. The colonoscopy, although reliable, is both relatively 
cosdy in money and time, and is an invasive, uncomfortable painful procedure for the 
patient 

20 The second detection technique is the xise of a barium enema and two- 

dimensional X-ray imaging of the colon. The barium enema is used to coat the colon 
with bariimi, and a two-dimensional X-ray image is taken to capture an image of the 
colon. However, barium enemas may not always provide a view of tiie entire colon, 
require extensive pretreatment and patient manipulation, is often operator-dependent 

25 ^en performing tibie operation, exposes the patient to excessive radiation and can be 
less sensitive than a colonoscopy. Due to deficiencies in the conventional practices 
described above, a more reliable, less intrusive and less expensive way to examine the 
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colon for polyps is desirable. A method to examine other human organs, such as the 
limgs, for masses in a reliable, cost effective way and with less patient discomfort is 
also desirable. 

Another leading cause of cancer deaths in the United States is bladder 

5 cancer. In 1995, there were 50,000 new cases of bladder cancer reported and 1 1,000 
deaths were reported as a result of this disease. The most common test for bladder 
cancer is the use of a urine "dipstick" or conventional urinalysis. However, such tests 
are generally only effective at detecting bladder cancer in its later developed stages 
and does not provide any information regarding the size or location of a cancerous 

1 0 growth. Cystoscopy, the mam method of investigating bladder abnormalities at 

present, provides accurate results and can provide information regarding the relative 
size and location of any abnormalities. However, cystoscopy is an invasive procedure 
which offers a physician a limited field of view and lacks an objective indication of 
size. In addition, cystoscopy is contra-indicated for those patients who have severe 

15 urethral strictures or active vesical bleeding. Thus, it is desirable to develop 

alternative procedures for screening patients for bladder cancer, especially at early 
stages of cancer development 

Two-dimensional ("2D") visualization of human organs employing 
currently available medical imaging devices, such as computed tomography and MRI 

20 (magnetic resonance imaging), has been widely used for patient diagnosis. Three- 
dimensional images can be formed by stacking and interpolating between two- 
dimensional pictures produced from the scanning machines. Imaging an organ and 
visualizing its volume in three-dimensional space would be beneficial due to its lack 
of physical intrusion and the ease of data manipulation. However, the exploration of 

25 the three-dimensional volume image must be properly performed in order to fiiUy 
e3q)loit the advantages of virtually viewing an organ from the inside. 

When viewing the three dimensional ("3D") volume virtual image of 
an enviroimient, a functional model must be used to explore the virtual space. One 
possible model is a vfrtual camera which can be used as a point of reference for the 

30 viewer to explore the virtual space. Camera control ia the context of navigation 

within a general 3D virtual envirormient has been previously studied. There are two 
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conventional types of camera control offered for navigation of virtual space. The jSrst 
gives the operator fiill control of the camera which allows the operator to manipulate 
the camera in different positions and orientations to achieve the view desired. The 
operator will in effect pilot the camera. This allows the operator to explore a 
5 particular section of interest while ignoring other sections. However, complete 
control of a camera in a large domain would be tedious and tiring, and an operator 
might not view all the important features between the start and finishing point of the 
exploration. 

The second technique of camera control is a planned navigation 
10 method, which assigns the camera a predetermined path to take and which cannot be 
changed by the operator. This is akin to having an engaged "autopilot". This allows 
the operator to concentrate on the virtual space being viewed, and not have to worry 
about steering into walls of the environment being examined. However, this second 
technique does not give the viewer the flexibility to alter the course or investigate an 
1 S interesting area viewed along the flight path. 

m 

It would be desirable to use a combination of the two navigation 
techniques described above to realize the advantages of both techniques while 
minimizing their respective drawbacks. It would be desirable to apply a flexible 
navigation technique to the examination of human or animal organs which are 

20 represented in virtual 3D space in order to perform a non-intrusive painless thorough 
examination. The desired navigation technique would further allow for a complete 
examination of a virtual organ in 3D space by an operator allowing flexibility while 
ensuring a smooth path and complete examination through and around the organ. It 
would be additionally desirable to be able to display the exploration of the organ in a 

25 real time setting by using a technique which minimizes the computations necessary 
for viewing the organ. The desired technique should also be equally applicable to 
exploring any virtual object. 

It is another object of the invention to assign opacity coefficients to 
each volume element in the representation in order to make particular volume 

30 elements transparent or translucent to varymg degrees in order to customize the 
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visualization of the portion of the object being viewed. A section of the object can 
also be composited using the opacity coefiBcients. 

SUMMARY OF THE INVENTION 
The invention generates a three-dimensional visualization image of an 
5 object such as a human organ using volume visualization techniques and explores the 
virtual image using a guided navigation system which allows the operator to travel 
along a predefined flight path and to adjust both the position and viewing angle to a 
particular portion of interest in the image away firom the predefined path in order to 
identify polyps, cysts or other abnormal features in the organ. 

10 A method for performing virtual examination of an object includes 

performing at least one imaging scan of an object with the object distended by the 
presence of a contrast agent. In addition, at least one imaging scan of the object is 
acquired with the object relieved of the contrast agent. The scans are converted to 
corresponding volume datasets formed with a plxu^ality of voxels. Image segmentation 

1 5 is then performed to classify the voxels of each scan into a plurality of categories. 
The volume datasets of each scan are registered to a common coordinate system. A 
displaying operation can then be performed where corresponding images at least two 
of the volume datasets are substantially simultaneously displayed. Virtual navigation 
operations performed in one of the volume datasets results in having the 

20 corresponding navigation operations take place in at least one other volume dataset. 

Preferably, the at least one scan of the distended object includes a 
transverse scan and a coronal scan of the object. Similarly, it is preferable that the at 
least one scan of the relieved object includes a transverse scan and a coronal scan of 
the object. This procedure is particularly well suited for performing virtual 

25 cystoscopy, where the object is the bladder. In this case, the scan generally takes the 
form of a magnetic resonance imaging scan and the contrast agent can be urine. 

Another method in accordance with the present invention is for 
performing virtual examination of an object. In this method, an imaging scan of the 
object is performed to acquire image scan data. The acquired image scan data is 

30 converted to a plurality of volume units, or voxels. By interpolating between the 
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voxels, an expanded dataset is generated. Image segmentation can then be performed 
to classify the voxels into a plurality of categories. A volume of the object interior is 
extracted from the expanded dataset, such as by using a region growing algorithm 
from a seed voxel within the object lumen. A reduced resolution dataset is then 
5 generated from the expanded dataset. To efficiently store and recall the data from the 
expanded data set, this data is stored in a tree data structure. Images can then be 
rendered for both the expanded dataset and reduced resolution dataset. One of these 
images is then selected for viewing. Generally, the reduced resolution dataset is 
selected during navigation or image interaction whereas the expanded dataset is 

10 selected for high resolution, static display. 

A method ofperforming virtual angiography is also provided. In this 
method, imaging scan data is acquired of at least a portion of the aorta. The imaging 
scan data is converted to a volume representation including a plurality of voxels. The 
volume representation is segmented to classify the voxels into one of a plurality of 

1 5 categories. The segmented volume representation is then analyzed to identify voxels 
indicative of at least a portion of an aneurysm in the aortic wall. From the portions of 
the aneurysm which are identified, at least one closing surface is generated around the 
voxels indicative of at least a portion of an aneurysm. The closing surface provides an 
estimate of the contour of the aneurysm. A navigation path can be established 

20 through the aortic lumen and characteristics of the aneurysm, such as length, diameter, 
volume and composition, can be determined. 

The method of performing virtual angiography can be used to detect 
and monitor the progression of aneurysms and can also be used in determining the 
characteristics needed to place a stent graft. 

23 Also provided is a method of defining a skeleton for a three 

dimensional image representation of a hollow object formed with a plurality of 
voxels. A root voxel is first identified within the hollow object. A distance map is 
then generated for all voxels within the hollow object. The distance map is formed 
using a 26*connected cubic plate having Euclidian weighted distances. Those voxels 

30 having a local maxima in the distance map are identified as endpoints of branches in 
the hollow object. For each local maxima voxel, a shortest connected path to either 
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the root voxel or a previously defined shortest path, is determined. The collection of 
shortest paths is the rough skeleton of the object This technique is particularly well 
suited for multibranch structures such as the respiratory system and cardio vascular 
system. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Further objects, features and advantages of the invention will become 
apparent from the following detailed description taken in conjunction with the 
accompanying figures showing a preferred embodiment of the invention, on which: 

Figure 1 is a flow chart of the steps for performing a virtual 
examination of an object, specifically a colon, in accordance with the invention; 

Figure 2 is an illustration of a "submarine" camera model which 
performs guided navigation in the virtual organ; 

Figure 3 is an illustration of a pendulum used to model pitch and roll of 
the "submarine" camera; 

Figure 4 is a diagram illustrating a two dimensional cross-section of a 
volumetric colon which identifies two blocking walls; 

Figure 5 is a diagram illustrating a two dimensional cross-section of a 
volumetric colon upon which start and finish volume elements are selected; 

Figure 6 is a diagram illustrating a two dimensional cross-section of a 
volumetric colon which shows a discrete sub-volume enclosed by the blocking walls 
and the colon surface; 

Figure 7 is a diagram illustrating a two dimensional cross-section of a 
volumetric colon which has multiple layers peeled away; 

Figure 8 is a diagram illustrating a two dimensional cross-section of a 
volumetric colon which contains the remaining flight path; 

Figure 9 is a flow chart of the steps of generating a volume 
visualization of the scanned organ; 

Figure 10 is an illustration of a vhtual colon which has been sub- 
divided into cells; 
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Figure 1 1 A is a graphical depiction of an organ which is being virtually 

examined; 

Figure 1 IB is a graphical depiction of a stab tree generated when 
depicting the organ in Fig. 1 1 A; 
5 Figure 1 1 C is a further graphical depiction of a stab tree generated 

while depicting the organ in Fig. 1 1 A; 

Figure 12A is a graphical depiction of a scene to be rendered with 
objects within certain cells of the scene; 

Figure 12B is a graphical depiction of a stab tree generated while 
1 0 depicting the scene in Fig. 12 A; 

Figures 12C-12E are further graphical depictions of stab trees 
generated while depicting the image in Fig. 1 2A; 

. Figure 1 3 is a two dimensional representation of a virtual colon 
containing a polyp whose layers can be removed; 
1 5 Figure 1 4 is a diagram of a system used to perform a virtual 

examination of a human organ in accordance with the invention; 

Figure 15 is a flow chart depicting an improved image segmentation 

method; 

Figure 16 is a graph of voxel intensity versus frequency of a typical 
20 abdominal CT data set; 

Figure 17 is a perspective view diagram of an intensity vector structure 
including a voxel of interest and its selected neighbors; 

Figure 1 8 A is an exemplary image slice from a CT scan of a human 
abdominal region, primarily illustrating a region including the lungs; 
25 Figure 18B is a pictorial diagram illustrating the identification of the 

lung region in the image slice of Figure 18 A; 

Figure 18C is a pictorial diagram illustrating the removal of the lung 
volume identified in Figure 18B; 

4 

Figure 19A is a exemplary image slice form a CT scan of a human 
30 abdominal region, primarily illustrating a region including a portion of the colon and 
bone; 



wo 02/29764 PCT/USOl/30704 

8 

Figure 19B is a pictorial diagram illustrating the identification of the 
colon and bone region from the image slice of Figure 1 9A; 

Figure 19C is a pictorial diagram illustrating the image scan of figure 
1 9a with the regions of bone removed; and 
5 Figure 20 is a flowchart illustrating a method for applying texture to 

monochrome image data. 

Figure 21 is a flowchart illustrating a method for volume rendering 

■ * 

employing a fast perspective ray castmg technique; 

Figure 22 is a flowchart illustrating a method for determining the 
10 central flight path through a colon lumen employing a volume shrinking technique. 

Figure 23 is a flowchart further illustrating a volume shrinking 
technique for use in the method illustrated in Figure 22. 

Figure 24 is a three dimensional pictorial representation of a 
segmented colon lumen with a central fly-path generated therein. 
15 - Figure 25 is a flow chart illustrating a method of generating a central 

flight path through a colon lumen employing a segmentation technique. 

Figure 26 is a block diagram of a system embodiment based on a 
personal computer bus architecture. 

Figure 27 is a flow chart illustrating a method of performing volume 

20 imaging using the system of Figure 26. 

Figure 28 is a flow chart illustrating a multi-scan method for 
performing virtual examination of an object, such as a bladder (virtual cystoscopy). 

Figure 29 is a pictorial representation of a display window suitable for 
presenting imaging results from the virtual cystoscopy method of Figure 28 and 
25 providing illustrative outside views of a bladder structure. • 

Figure 30 is a pictorial representation of a display window suitable for 
presenting imaging results from the virtual cystoscopy method of Figure 28 and 
providing illustrative interior views of a bladder structure. 

Figure 31 is a flow chart of a method of performing virtual 
30 examination of an object, such as the larynx, using multiresolution viewing. 
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Figure 32 is a flow chart of a method for performing virtual 

angiography. 

Figures 33 A-C are pictorial views of a portion of the aorta illustrating 
the presence of an abdominal aortic aneurysm. 
5 Figure 34 is a flow chart illustrating a method for generating a skeleton 

structure of an object. 

Figure 35 is a schematic diagram of a 26-coimected, Euclidean 
weighted, cubic distance plate. 

Figure 36 is a diagram illustrating pseudo-code of a process for 
0 generating a distance map for use in the method of Figure 34. 

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 

While the methods and systems described in this application can be 
applied to any object to be examined, the preferred embodiment which will be 
described is the examination of an organ in the human body, specifically the colon. 
5 The colon is long and twisted which makes it especially suited for a virtual 
examination saving the patient both money and the discomfort and danger of a 
physical probe. Other examples of organs which can be examined, without limitation, 
include the lungs, stomach and portions of the gastro-intestinal system, the heart and 
blood vessels. 

0 Fig. 1 illustrates the steps necessary to perform a virtual colonoscopy 

using volume visualization techniques. Step 101 prepares the colon to be scanned in 
order to be viewed for examination if required by either the doctor or the particular 
scanning instrument. This preparation could include cleansing the colon with a 
"cocktail" or liquid which enters the colon after being orally ingested and passed 

5 through the stomach. The cocktail forces the patient to expel waste material that is 
present in the colon. One example of a substance used is Golytely. Additionally, in 
the case of the colon, air or CO2 can be forced into the colon in order to expand it to 
make the colon easier to scan and examine. This is accomplished with a small tube 
placed in the rectum with approximately 1 ,000 cc of air pumped into the colon to 

iO distend the colon. Depending upon the type of scanner used, it may be necessary for 
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the patient to drink a contrast substance such as barium to coat any unexpunged stool 
in order to distinguish the waste in the colon from the colon walls themselves. 
Alternatively, the method for virtually examining the colon can remove the virtual 
waste prior to or during the virtual examination as explained later in this specification. 
5 Step 101 does not need to be performed in all examinations as indicated by the dashed 
line in Fig. 1 . 

Step 103 scans the organ which is to be examined. The scanner can be 
an apparatus well known in the art, such as a spiral CT-scanner for scanning a colon 
or a Zenita MRI machine for scanning a lung labeled for example with xenon gas. 

0 The scanner must be able to take multiple images from different positions around the 
body during suspended respiration, in order to produce the data necessary for the 
volume visualization. An example of a single CT-image would use an X-ray beam of 
5mm width, 1 :1 to 2:1 pitch, with a 40cm field-of-view being performed from the top 
of the splenic flexure of the colon to the rectum. 

S Discrete data representations of said object can be produced by other 

methods besides scanning. Voxel data representing an object can be derived from a 
geometric model by techniques described in U.S. Pat. No. 5,038,302 entitled "Method 
of Converting Continuous Three-Dimensional Geometrical Representations into 
Discrete Three-Dimensional Voxel-Based Representations Within a Three- 

0 Dimensional Voxel-Based System" by Kaufman, issued Aug. 8, 1991, filed July 26, 
1 988, which is hereby incorporated by reference. Additionally, data can be produced 
by a computer model of an image which can be converted to three-dimension voxels 
and explored in accordance with this invention. One example of this type of data is a 
computer simulation of the turbulence surrounding a space shuttle craft. 

.5 Step 1 04 converts the scanned images into three-dimensional volume 

elements (Voxels). Li the preferred embodiment for examining a colon, the scan data 
is reformatted into 5mm thick slices at increments of 1mm or 2.5mm and 
reconstructed in 1mm slices, with each slice represented as a matrix of 5 12 by 
512 pixels. By doing this, voxels of approximately 1 cubic mm are created. Thus a 

0 large nimiber of 2D slices are generated depending upon the length of the scan. The 
set of 2D slices is then reconstructed to 3D voxels. The conversion process of 2D 
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images from the seamier into 3D voxels can either be performed by the scanning 
machine itself or by a separate machine such as a computer with techniques which are 
well known in the art (for example, see U.S. Pat. No. 4,985,856 entitled "Method and 
" Apparatus for Storing, Accessing, and Processing Voxel-based Data" by Kaufinan et 
5 al.; issued Jan. 15, 1991, filed Nov. 11, 1988; which is hereby incorporated by 
reference). 

Step 1 05 allows the operator to define the portion of the selected organ 
to be examined. A physician may be interested in a particular section of the colon 
likely to develop polyps. The physician can view a two dimensional slice overview 

10 map to indicate the section to be examined. A starting point and finishing point of a 
path to be viewed can be indicated by the physician/operator. A conventional 
computer and computer interface (e.g., keyboard, mouse or spaceball) can be used to 
designate the portion of the colon which is to be inspected. A grid system with, 
coordinates can be used for keyboard entry or the physician/operator can "click" on 

15 the desired points. The entire image of the colon can also be viewed if desired. 

Step 1 07 performs the planned or guided navigation operation of the 
virtual organ being exammed. Performing a guided navigation operation is defined as 
navigating through an environment along a predefined or automatically predetermined 
flight path which can be manually adjusted by an operator at any time. After the scan 

20 data has been converted to 3D voxels, the inside of the organ must be traversed from 
the selected start to the selected finishing point. The virtual examinations is modeled 
on having a tiny camera traveling through the virtual space with a lens pointing 
towards the finishing point The guided navigation technique provides a level of 
interaction with the camera, so that the camera can navigate through a virtual 

25 environment automatically in the case of no operator interaction, and at the same time, 
allow the operator to manipulate the camera when necessary. The preferred 
embodiment of achieving guided navigation is to use a physically based camera model 
which employs potential fields to control the movement of the camera and which are 
described in detail in Figs. 2 and 3. 

30 Step 109, which can be performed concurrently with step 107, displays 

the inside of the organ from the viewpoint of the camera model along the selected 
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pathway of the guided navigation operation. Three-dimensional displays can be 
generated using techniques well known in the art such as the marching cubes . 
technique. However, in order to produce a real time display of the colon, a technique 
is required which reduces the vast number of computations of data necessary for the 
5 display of the virtual organ. Fig. 9 describe this display step in more detail. 

The method described in Figure 1 can also be applied to scanning 
multiple organs in a body at the same time. For example, a patient may be examined 
for cancerous growths in both the colon and lungs. The method of Figure 1 would be 
modified to scan ail the areas of interest in step 103 and to select the current organ to 

10 be examined in step 105. For example, the physician/operator may initially select the 
colon to virtually explore and later explore the lung. Alternatively, two different 
doctors with different specialties may virtually explore different scanned organs 
relating to their respective specialties. Following step 109, the next organ to be 
examined is selected and its portion will be defined and explored. This continues 

IS until all organs which need examination have been processed. 

The steps described in conjunction with Figure 1 can also be applied to 
the exploration of any object which can be represented by volume elements. For 
example, an architectural structure or inanimate object can be represented and 
explored in the same manner. 

20 Figure 2 depicts a "submarine" camera control model which performs 

the guided navigation technique in step 107. When there is no operator control during 
guided navigation, the default navigation is similar to that of planned navigation 
which automatically directs the camera along a flight path fi'om one selected end of 
the colon to another. During the plarmed navigation phase, the camera stays at the 

25 center of the colon for obtaining better views of the colonic surface. When an 
interesting region is encountered, the operator of the virtual camera using guided 
navigation can interactively bring the camera close to a specific region and direct the 
motion and angle of the camera to study the interesting area in detail, without 
imwillingly colliding with the walls of the colon. The operator can control the camera 

30 with a standard interface device such as a keyboard, mouse or non-standard device 
such as a spaceball. In order to fully operate a camera in a virtual environment, six 
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degrees of freedom for the camera is required. The camera must be able to move in 
the horizontal, vertical, and Z direction (axes 217), as well as being able to rotate in 
another three degrees of freedom (axes 219) to allow the camera to move and scan all 
sides and angles of a virtual environment. The camera model for guided navigation 
5 includes an inextensible, weightless rod 201 connectmg two particles x, 203 and Xj 
205, both particles being subjected to a potential field 215. The potential field is 
defined to be highest at the walls of the organ in order to push the camera away from 
the walls. 

The positions of the particles are given by x, and Xj, and they are 
[) assumed to have the same mass m. A camera is attached at the head of the submarine 
X, 203, whose viewing direction coincides with XjX,. The submarine can perform 
translation and rotation around the center of mass x of the model as the two particles 
are affected by the forces from the potential field V(x) which is defined below, any 
friction forces, and any simulated extemal force. The relations between x„ Xj, and x 
5 are as follows: 

. r = (rsinOcos%rsinOsin%rcosO)y 
X, = X + r, 

Xj = X - r, (1) 



I) 



where r, 6 and 9 are the polar coordinates of the vector xXj. 

The kinetic energy of the model, T, is defined as the summation of the kinetic 

energies of the movements of x, and X2: 

= mx + mr 

= rn(i^ + y^ + i') + 7nr^(^^+^2sin^<?). (2) 

Then, the equations for the motion of the submarine model are 
S obtained by using LaGrange's equation: 
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(3) 



where the qjS are the generalized coordinates of the model and can be considered as 
the variables of time / as: 



with v denoting the roll angle of our camera system, which will be explained later. 
. The FjS are called the generalized forces. The control of the submarine is performed 
5 by applying a simulated external force to x„ 



and it is assumed that both Xi and x^ are affected by the forces from the potential field 
and the frictions which act in the opposite direction of each particle's velocity. Conse 
quently, the generalized forces are formulated as follows: 



10 where k denotes the friction coefficient of the system. The external force Fe,, is 

applied by the operator by simply clicking the mouse button in the desired direction 
207 in the generated image, as shown in Figure 2. This camera model would then be 
moved in that direction. This allows the operator to control at least five degrees of 
freedom of the camera with only a single click of the mouse button. From Equations 

15 (2), (3) and (5), it can be derived that the accelerations of the five parameters of our 
submarine model as: 



(4) 




Pi = -mVV(xi)-Axi + ' esti 



(5) 
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= _l/ gV^(xi) , aV(x,) , fci . 



XT = 

9 = ^'dji0cos9 



-(cos tf{co*^(— ^ - + sm^(-5--i - -^)} 



2r' 02 5« ' ^ 5j/ 

- m 2mr* 

sin 5- ) + coa^(— r- S— )> 



2r^ ^> fix 5x ' ^ 

. 1 

^ ^ (-f^ sin ^ + cos ^)], (6) 



m 2mr 



where x and x denote the first and the second derivative of x, respectively, and 



j' avU)_ ^ dV(x) ^ dV(x) \ jignotes the gradient of the potential at a point x. 
^ 8x dy az J 



• ■ 

The terms (p^sinecos9 of 9 and _^£$_£££_ of $ are called the 

sin6 



centrifugal force and the Coriolis force^ respectively, and they are concerned with the 
5 exchange of singular velocities of the submarine. Since the model does not have the 
moment of inertia defined for the rod of the submarine, these terms tend to cause an 
overflow of the nimieric calculation of <p. Fortunately, these terms become significant 
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only when the angular velocities of the submarine model are significant, which 
essentially means that the camera moves too fast. Since it is meaningless to allow the 
camera to move so fast because the organ could not be properly viewed, these terms 
are minimized in our implementation to avoid the overflow problem. 
5 From the first three formulas of Equation (6), it is known that the 

submarine cannot be propelled by the external force against the potential field if the 
following condition is satisfied: 

lF„t| 



|VV(xx) + VV(xOI> 



m 



Since the velocity of the submarine and the external force have upper limits in our 

implementation, by assigning sufficiently high potential values at the boundary of the 
10 objects, it can be guaranteed that the submarine never bumps against the objects or 

walls in the environment. 

As mentioned previously, the roll angle v|/ of the camera system needs 

to be considered. One possible option allows the operator fiill control of the angle \|;. 

However, although the operator can rotate the camera fi-eely around the rod of the 
1 S model, he or she can easily become disoriented. The preferred technique assumes that 

the upper direction of the camera is connected to a pendulum with mass 301, 

which rotates fi-eely around the rod of the submarine, as shown in Figure 3. The 

direction of the pendulum, rj, is expressed as: 



r2 = r2(cos tfcos^sinV' + sin <^ cos V^, cos 9 sin ^sin — cos ^cos V^i sin ^ sin •0). 



although it is possible to calculate the accurate movement of this pendulum along 
20 with the movement of the submarine, it makes the system equations too complicated. 
Therefore, it is assumed that all the generalized coordinates except the roll angle \|/ are 
constants, and thus define the independent kinetic energy for the pendulum system as: 
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7112 2 mmnL-'h^ 

This simplifies the model for the roll angle. Since it is assumed in this model that the 
gravitational force 

= m2g = (m2^x>7n2^y,m25x) 

acts at the mass point the acceleration of \|f can be derived using LaGrange*s 
equation as: 

- 1 

^ = — {{^aj(cos5cos^cos^ — sin^sinV^) 
+gy{cos0 wi^cosij} + cos^sinV^) 

siB^JcosV^)} ^. (7) 

77l2 

5 From Equations (6) and (7), the generalized 

coordinates q(t) and their derivatives q(t) are calculated asymptotically by using 
Taylor series as: 

Hit + h) = q{t) + hii{t) + jq{t) + 0{h% 
q{t + k) = q{t) + hqit) + 0{h% 

to freely move the submarine. To smooth the submarine's motion, the time step h is 
10 selected as an equilibrium value between being as small as possible to smooth the 
motion but as large as necessary to reduce computation cost. 



I 
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Definition of the Potential Field 
The potential field in the submarine model in Figure 2 defines the 



10 



15 



20 



25 



boundaries (walls or other matter) in the virtual organ by assigning a high potential to 
the boundary in order to ensure that the submarine camera does not collide with the 
walls or other boundary. If the camera model is attempted to be moved into a high 
potential area by the operator, the camera model will be restrained from doing so 
unless the operator wishes to examine the organ behind the boundary or inside a 
polyp, for example. In the case of performing a virtual colonoscopy, a potential field 
value is assigned to each piece of volxmietric colon data (volume element). When a 
particular region of interest is designated in step 105 of Fig. 1 with a start and finish 
point, the voxels within the selected area of the scanned colon are identified using 
conventional blocking operations. Subsequently, a potential value is assigned to every 
voxel x of the selected volume based on the following three distance values: the 
distance firom the finishing point dt(x), the distance firom the colon surface ds(x) and 
the distance firom the center-line of the colon space dc(x). dt(x) is calculated by using 
a conventional growing strategy. The distance from the colon surface, ds(x), is 
computed using a conventional technique of growing fi:om the surface voxels inwards. 
To determine dc(x), the center-line of the colon fi-om the voxel is first extracted, and 
then dc(x) is computed using the conventional growing strategy from the center-line 
of the colon. 

To calculate the center-line of the selected colon area defined by the 
user-specified start point and the user-specified finish point, the maximum value of 
ds(x) is located and denoted dmax. Then for each voxel inside the area of interest, a 
cost value of dmax - ds(x) is assigned. Thus the voxels which are close to the colon 
surface have high cost values and the voxels close to the center line have relatively 
low cost values. Then, based on the cost assignment, the single-source shortest path 
technique which is well known in the art is applied to efficiently compute a minimum 
cost path from the source point to the finish point. This low cost line indicates the 
center-line or skeleton of the colon section which is desired to be explored. This 
technique for determining the center-line is the preferred technique of the invention. 



# 
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To compute the potential value V(x) for a voxel x inside the area of 
interest, the following formula is employed: 

vw = c,*(xr + c,(;^^i^r, (8) 

where C|, Cj, n and v are constants chosen for the task. In order to avoid any 
collision between the virtual camera and the virtual colonic surface, a sufficiently 
5 large potential value is assigned for all points outside the colon. The gradient of the 
potential field will therefore become so significant that the submarine model camera 
will never collide with the colonic wall when being run. 

Another technique to determine the center-line of the path in the colon 
is called the "peel-layer" technique and is shown in Figure 4 through Figure 8. 

10 Figure 4 shows a 2D cross-section of the volumetric colon, with the 

two side walls 401 and 403 of the colon being shown. Two blocking walls are 
selected by the operator in order to define the section of the colon which is of interest 
to examine. Nothing can be viewed beyond the blocking w^ls. This helps reduce the 
number of computations when displaying the virtual representation. The blocking 

1 5 walls together ^yith side walls identify a contained voliunetric shape of the colon 
which is to be explored. 

Figure 5 shows two end points of the flight path of the virtual 
examination, the start volume element 501 and the fmish volume element 503. The 
start and finish points are selected by the operator in step 1 05 of Fig. 1 . The voxels 

20 between the start and finish points and the colon sides are identified and marked, as 
indicated by the area designated with "x"s in Fig. 6. The voxels are three-dimensional 
representations of the picture element 

The peel-layer technique is then applied to the identified and marked 
voxels in Fig. 6. The outermost layer of all the voxels (closest to the colon walls) is 

25 peeled off step-by-step, until there is only one inner layer of voxels remaining. Stated 
differently, each voxel furthest away firom a center point is removed if the removal 
does not lead to a disconnection of the path between the start voxel and the finish 
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voxel. Figure 7 shows the intermediate result after a number of iterations of peeling 
the voxels in the virtual colon are complete. The voxels closest to the walls of the 
colon have been removed. Fig. 8 shows the final flight path for the camera model 
down the center of the colon after all the peeling iterations are complete. This 
5 produces essentially a skeleton at the center of the colon and becomes the desired 
flight path for the camera model. 

Buffer assisted visibiUtv 
Figure 9 describes a real time visibility technique to display of virtual 
images seen by the camera model in the virtual three-dimensional volume 

0 representation of an organ. Figure 9 shows a display technique using a modified Z 
buffer which corresponds to step 109 in Fig. 1. The number of voxels which could be 
possibly viewed from the camera model is extremely large. Unless the total number 
of elements (or polygons) which must be computed and visualized is reduced fi'om an 
entire set of voxels in the scanned environment, the overall number of computations 

5 will make the visualization display process exceedingly slow for a large internal area. 
However, in the present invention only those images which are visible on the colon 
surface need to be computed for display. The scanned environment can be subdivided 
into smaller sections, or cells. The Z buffer technique then renders only a portion of 
the cells which are visible fi-om the camera. The Z buffer technique is also used for 

0 three-dimensional voxel representations. The use of a modified Z buffer reduces the 
number of visible voxels to be computed and allows for the real time examination of 
the virtual colon by a physician or medical technician. 

The area of interest firom which the center-line has been calculated in 
step 107 is subdivided into cells before the display technique is applied. Cells are 

:5 collective groups of voxels which become a visibility unit. The voxels in each cell 
will be displayed as a group. Each cell contains a number of portals through which 
the other cells can be viewed. The colon is subdivided by beginning at the selected 
start point and movmg along the center-line 1001 towards the finish point. The colon 
is then partitioned into cells (for example, cells 1003, 1005 and 1007 in Fig. 10) when 

10 a predefined threshold distance along the center-path is reached. The threshold 
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distance is based upon the specifications of the platform upon which the visuaUzation 
technique is perfonned and its capabilities of storage and processing. The cell size is 
directly related to the number of voxels which can be stored and processed by the 
platfonn. One example of a threshold distance is 5cm, although the distance can 

5 greatly vary. Each cell has two cross-sections as portals for viewing outside of the 
cell as shown in Fig. 10. 

Step 901 in Fig. 9 identifies the cell within the selected organ which 
currently contains the camera. The current cell will be displayed as well as all other 
cells which are visible given the orientation of the camera. Step 903 builds a stab tree 

10 (tree diagram) of hierarchical data of potentially visible cells from the camera 

(through defined portals), as will be described in fiirther detail hereinbelow. The stab 
tree contains a node for every cell which may be visible to the camera. Some of the 
cells may be transparent without any blocking bodies present so that more than one 
cell will be visible in a single direction. Step 905 stores a subset of the voxels firom a 

15 cell which include the intersection of adjoining cell edges and stores them at the 
outside edge of the stab tree in order to more efficiently determine which cells are 
visible. 

Step 907 checks if any loop nodes are present in the stab tree. A loop 
node occurs when two or more edges of a single cell both border on the same nearby 
20 cell. This may occur when a single cell is surrounded by another cell. If a loop node 
is identified in the stab tree, the method continues with step 909. If there is no loop 
node, the process goes to step 911. 

Step 909 collapses the two cells making up the loop node into one 
large node. The stab tree is then corrected accordingly. This eliminates the problem 
25 of viewing the same cell twice because of a loop node. The step is performed on all 
identified loop nodes. The process then continues with step 911. 

Step 91 1 then initiates the Z-buffer with the largest Z value. The Z 
value defines the distance away from the camera along the skeleton path. The tree is 
then traversed to first check the intersection values at each node. If a node 
30 intersection is covered, meaning that the current portal sequence is occluded (which is 
determined by the Z buffer test), then the traversal of the current branch in the tree is 



wo 02/29764 PCT/USOl/30704 

22 

stopped. Step 913 traverses each of the branches to check if the nodes are covered 
and displays them if they are not. 

Step 915 then constructs the image to be displayed on the operator's 
screen from the volume elements within the visible cells identified in step 913 using 

) one of a variety of techniques known in the art, such as volume rendering by 
compositing. The only cells shown are those which are identified as potentially 
visible. This technique limits the number of cells which requires calculations in order 
to achieve a real time display and correspondingly increases the speed of the display 
for better performance. This technique is an improvement over prior techniques 

) which calculate all the possible visible data points whether or not they are actually 
viewed. 

Figure 1 1 A is a two dimensional pictorial representation of an organ 
which is being explored by guided navigation and needs to be displayed to an 
operator. Organ 1 101 shows two side walls 1 102 and an object 1 105 in the center of 

5 the pathway. The organ has been divided into four cells A1151,B1153, C1155 and 
D 1 157. The camera 1 103 is facing towards cell D 1 157 and has a field of vision 
defined by vision vectors 1 107, 1 108 which can identify a cone-shaped field. The 
cells which can be potentially viewed are cells B 1153, C 1 155 and D 1 157. Cell C 
1 155 is completely surrounded by Cell B and thus constitutes a node loop. 

[) Fig. 1 IB is a representation of a stab tree built firpm the cells in Fig. 

1 1 A. Node A 1 109 which contains the camera is at the root of the tree. A sight line 
or sight cone, which is a visible path without being blocked, is drawn to node B 1110. 
Node B has direct visible sight lines to both node C 1 1 12 and node D 1 1 14 and which 
is shown by the connecting arrows. The sight line of node C 1 1 12 in the direction of 

5 the viewing camera combines with node B 1 1 1 0. Node C 1 1 12 and node B 1 1 10 will 
thus be collapsed into one large node B' 1 122 as shown in Fig. 1 IC. 

Fig. 1 IC shows node A 1 109 containing the camera adjacent to node B* 
1 122 (containing both nodes B and node C) and node D 1 1 14. The nodes A, B' and D 
will be displayed at least partially to the operator. 

0 Figs 12A - 12E illustrate the use of the modified Z buffer with cells 

that contain objects which obstruct the views. An object could be some waste 
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material in a portion of the virtual colon. Fig. 12A shows a virtual space with 
10 potential cells: A 1251, B 1253, C 1255, D 1257, E 1259, F 1261. G 1263, H 1265, 
1 1267 and J 1269. Some of the cells contain objects. If the camera 1201 is 
positioned in cell 1 1267 and is facing toward cell F 1261 as indicated by the vision 

5 vectors 1203, then a stab tree is generated in accordance with the technique illustrated 
. by the flow diagram in Fig. 9. Fig. 12B shows the stab tree generated with the 
intersection nodes showing for the virtual representation as shown in Fig. 12 A. Fig. 
12B shows cell 1 1267 as the root node of the tree because it contams the camera 
1201. Node 1 121 1 is pointing to node F 1213 (as indicated with an arrow), because 

10 cell F is directly connected to the sight line of the camera. Node F 1213 is pointing to 
both node B 1215 and node E 1219. Node B 1215 is pointing to node A 1217. Node 
C 1202 is completely blocked from the line of sight by camera 1201 so is not included 
in the stab tree. 

Fig. I2C shows the stab tree after node 1 121 1 is rendered on the 
1 5 display for the operator. Node 1 1 2 1 1 is then removed from the stab tree because it 
has abready been displayed and node F 121 3 becomes the root. Fig. 12D shows that 
node F 12 1 3 is now rendered to join node 1.1211. The next nodes in the tree 
connected by arrows are then checked to see if they are already covered (already 
processed). In this example, all of the intersected nodes from the camera positioned 
20 in cell 1 1267 has been covered so that node B 5 15 (and therefore dependent node A) 
do not need to be rendered on the display. 

Fig. 12E shows node E 515 being checked to determine if its 
intersection has been covered. Since it has, the only rendered nodes in this example 
of Figure 12A-1 2E are nodes I and F while nodes A, B and E are not visible and do 
25 not need to have their cells prepared to be displayed. 

The modified Z buffer technique described in Figure 9 allows for fewer 
computations and can be applied to an object which has been represented by voxels or 
other data elements, such as polygons. 

Figure 13 shows a two dimensional virtual view of a colon with a large 
30 polyp present along one of its walls. Figure 1 3 shows a selected section of a patient's 
colon which is to be examined further. The view shows two colon walls 1301 and 
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1303 with the growth indicated as 1305. Layers 1307, 1309, and 131 1 show inner 
layers of the growth. It is desirable for a physician to be able to peel the layers of the 
polyp or tumor away to look inside of the mass for any cancerous or other harmful 
material. This process would in effect perform a vutual biopsy of the mass without 
5 actually cutting into the mass. Once the colon is represented virtually by voxels, the 
process of peeling away layers of an object is easily performed in a similar manner as 
described in conjunction with Figs. 4 through 8. The mass can also be sliced so that a 
particular cross-section can be examined. In Fig. 1 3 , a planar cut 1 3 1 3 can be made 
so that a particular portion of the growth can be examined. Additionally, a user- 
1 0 defined slice 1319 can be made in any manner in the growth. The voxels 1 3 1 9 can 

■ 

either be peeled away or modified as explained below. 

A transfer function can be performed to each voxel in the area of 
interest which can make the object transparent, semi-transparent or opaque by altering 
coefficients representing the translucently for each voxel. An opacity coefficient is 

1 5 assigned to each voxel based on its density. A mapping function then transforms the 
density value to a coefficient representing its translucency. A high density scaimed 
voxel will indicate either a wall or other dense matter besides simply open space. An 
operator or program routine could then change the opacity coefficient of a voxel or 
group of voxels to make them appear transparent or semi-transparent to the submarine 

20 camera model. For example, an operator may view a tumor within or outside of an 
entire growth. Or a transparent voxel will be made to appear as if it is not present for 
the display step of Figure 9. A composite of a section of the object can be created 
using a weighted average of the opacity coefficients of the voxels in that section. 

If a physician desires to view the various layers of a polyp to look for a 

25 cancerous areas, this can be performed by removing the outer layer of polyp 1305 . 
yielding a first layer 1307. Additionally, the first inner layer 1307 can be stripped 
back to view second inner layer 1309. The second inner layer can be stripped back to 
view third inner layer 131 1, etc. The physician could also slice the polyp 1305 and 
view only those voxels within a desired section. The slicing area can be completely 

30 user-defined. 



wo 02/29764 PCT/USOl/30704 

■ 

25 

Adding an opacity coefficient can also be used in other ways to aid in 
the exploration of a virtual system. If waste material is present and has a density as 
other properties within a certain known range, the waste can be made transparent to 
the virtual camera by changing its opacity coefficient during the examination. This 
5 will allow the patient to avoid ingesting a bowel cleansing agent before ttie procedure 
and make the examination faster and easier. Other objects can be similarly made to 
disappear depending upon the actual appUcation. Additionally, some objects like 

* 

polyps could be enhanced electronically by a contrast agent followed by a use of an 
appropriate transfer function. * 

1 0 Figure 1 4 shows a system for performing the virtual examination of an 

object such as a human organ using the techniques described in this specification. 
Patient 1401 Ues down on a platform 1402 while scanning device 1405 scans the area 
that contains the organ or organs which are to be examined. The scanning device 
1405 contains a scanning portion 1403 which actually takes images of the patient and 

15 an electronics portion 1406. Electronics portion 1406 comprises an interface 1407, a 
central processing unit 1409, a memory 141 1 for temporarily storing the scanning 
data, and a second interface 1413 for sending data to the virtual navigation platform. 
Interface 1407 and 1413 could be included in a single interface component or could be 
the same component. The components in portion 1406 are connected together with 

20 conventional connectors. 

In system 1400, the data provided from the scanning portion of device 
1403 is transferred to portion 1405 for processing and is stored in memoiy 141 1. 
Central processing imit 1 409 converts the scanned 2D data to 3D voxel data and 
stores the results in another portion of memory 1411. Alternatively, the converted 

25 data could be directly sent to interface unit 1413 to be transferred to the virtual 

navigation terminal 1416. The conversion of the 2D data could also take place at the 
virtual navigation terminal 1416 after being transmitted from interface 1413. In the 
preferred embodiment, the converted data is transmitted over carrier 1414 to the 
virtual navigation terminal 1416 in order for an operator to perform the virhial 

30 examination. The data could also be transported in other conventional ways such as 
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Storing the data on a storage medium and physically transporting it to terminal 1416 or 
by using satellite transmissions. 

The scanned data may not be converted to its 3D representation until 
the visualization rendering engine requires it to be in 3D form. This saves 
5 computational steps and memory storage space. 

Virtual navigation temiinal 1416 includes a screen for viewing the 
virtual organ or other scanned image, an electronics portion 1415 and interface control 
1419 such as a keyboard, mouse or spaceball. Electronics portion 1415 comprises a 
interface port 1421, a central processing unit 1423, other components 1427 necessary 

10 to run the terminal and a memory 1425. The components in terminal 1416 are 
connected together with conventional connectors. The converted voxel data is 
received in interface port 1421 and stored in memory 1425. The central processor 
unit 1423 then assembles the 3D voxels into a virtual representation and runs the 
submarine camera model as described in Figures 2 and 3 to perform the virtual 

1 5 examination. As the submarine camera travels through the virtual organ, the visibility 
technique as described in Fig[ure 9 is used to compute only those areas which are 
visible from the virtual camera and displays them on screen 141 7. A graphics ' 
accelerator can also be used in generating the representations. The operator can use 
interface device 1419 to indicate which portion of the scanned body is desired to be 

20 explored. The interface device 1419 can further be used to control and move the 
submarine camera as desired as discussed in Figure 2 and its accompanying 
description. Terminal portion 1415 can be the Cube-4 dedicated system box, 
generally available from the Department of Computer Science at the State University 
of New York at Stony Brook. 

25 Scanning device 1405 and terminal 1416, or parts thereof, can be part 

of the same imit. A single platform would be used to receive the scan image data, 
connect it to 3D voxels if necessary and perform the guided navigation. 

An important feature in system 1400 is that the virtual organ can be 
examined at a later time without the presence of the patient Additionally, the virtual 

30 examination could take place while the patient is being scanned. The scan data can 
also be sent to multiple terminals which would allow more than one doctor to view 
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the inside of the organ simultaneously. Thus a doctor in New York could be looking 
at the same portion of a patient's organ at the same time with a doctor in California 
while discussing the case. Alternatively, the data can be viewed at different times. 
Two or more doctors could perform their own examination of the same data in a 
5 diflBcult case. Multiple virtual navigation terminals could be used to view the same 
scan data. By reproducing the organ as a virtual organ with a discrete set of data, 
there are a multitude of benefits in areas such as accuracy, cost and possible data 
manipulations. 

The above described techniques can be further enhanced in virtual 

10 colonoscopy applications through the use of an improved electronic colon cleansing 
technique which employs modified bowel preparation operations followed by image 
segmentation operations, such that fluid and stool remaining in the colon during a 
computed tomographic (CT) or magnetic resonance imaging (MRI) scan can be 
detected and removed from the virtual colonoscopy unages. Through the use of such 

1 5 techniques, conventional physical washing of the colon, and its associated 
inconvenience and discomfort, is minimized or completely avoided. 

Referring to Figure 15, the first step in electronic colon cleansing is 
bowel preparation (step 1510), which takes place prior to conducting the CT or 
magnetic resonance imaging (MRI) scan and is intended to create a condition where 

20 residual stool and fluid remaining in the colon present significantly different image 
properties from that of the gasrfiUed colon interior and colon wall. An exemplary 
bowel preparation operation includes ingesting three 250 cc doses of Barium Sulfate 
suspension of 2.1 % WA/^, such as manufactured by E-Z-EM, Inc.,of Westbury, New 
York, during the day prior the CT or MRI scan. The three doses should be spread out 

25 over the course of the day and can be ingested along with three meals, respectively. 
The Barium Sulfate serves to enhance the images of any stool which remains in the 
colon. In addition to the intake of Barium Sulfate, fluid intake is preferably increased 
during the day prior to the CT or MRI scan. Cranberry juice is known to provide 
increased bowel fluids and is preferred, although water can also be ingested. In both 

30 the evening prior to the CT scan and the morning of the CT scan, 60 ml of a 

Diatrizoate Meglumine and Diaztrizoate Sodium Solution, which is commercially 
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available as MD-Gastroview, manufactured by Mallinckrodt, Inc. of St. Louis, 
Missouri, can be consumed to enhance image properties of the colonic fluid. Sodium 
phosphate can also be added to the solution to liquidize the stool in the colon, which 
provides for more uniform enhancement of the colonic fluid and residual stool. 
5 The above described exemplary preliminary bowel preparation 

operation can obviate the need for conventional colonic washing protocols, which can 
call for the ingestion of a gallon of Golytely solution prior to a CT scan. 

Just prior to conductmg the CT scan, an intravenous injection of 1 ml 
of Glucagon, manufactured by Ely Lily and Company, of Indianapolis, Indiana can be 

10 administered to minimize colon collapse. Then, the colon can be inflated using 
approximately lOOOcc of compressed gas, such as COj, or room air, which can be 
introduced through a rectum tube. At this point, a conventional CT scan is performed 
to acquire data from the region of the colon (step 1 520). For example, data can be 
acquired using a GE/CTI spiral mode scanner operating in a helical mode of 5mm, 

15 1.5-2.0:1 pitch, reconstructed in 1mm slices, where the pitch is adjusted based upon 

« 

the patient's height in a known manner. A routine imaging protocol of 120 kVp and 
200-280 ma can be utilized for this operation. The data can be acquired and 
reconstructed as Imm thick slice images having an array size of 512x512 pixels in the 
field of view, which varies from 34 to 40 cm depending on the patient's size, the 

20 number of such slices generally varies under these conditions from 300 to 450, 
depending on the patient's height. The image data set is converted to volume 
elements or voxels (step 1 530). 

Image segmentation can be performed in a number of ways. In one 
present method of image segmentation, a local neighbor technique is used to classify 

25 voxels of the image data in accordance with similar intensity values. In this method, 
each voxel of an acquired image is evaluated with respect to a group of neighbor 
voxels. The voxel of interest is referred to as the central voxel and has an associated 
intensity value. A classification indicator for each voxel is established by comparing 
the value of the central voxel to each of its neighbors. If the neighbor has the same 

30 value as the central voxel, the value of the classification indicator is incremented. 
However, if the neighbor has a different value from the central voxel, the 
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classification indicator for the central voxel is decremented. The central voxel is then 
classified to that category which has the maximum indicator value, which indicates 
the most uniform neighborhood among the local neighbors. Each classification is 
indicative of a particular intensity range, which in turn is representative of one or 
5 more material types being imaged. The method can be further enhanced by 

employing a mixture probability function to the similarity classifications derived. 

An alternate process of image segmentation is performed as two major 
operations: low level processing and high level feature extraction. During low level 
processing, regions outside the body contour are eliminated from further processing 

1 0 and voxels within the body contour are roughly categorized in accordance witti well 
defined classes of intensity characteristics. For example, a CT scan of the abdominal 
region generates a data set which tends to exhibit a well defined intensity distribution. 
The graph of Figure 1 6 illustrates such an intensity distribution as an exemplary 
histogram having four, well defined peaks, 1602, 1604, 1606, 1608, which can be 

1 S classified according to intensity thresholds. 

The voxels of the abdominal CT data set are roughly classified as four 
clusters by intensity thresholds (step 1540). For example, Cluster 1 can include 
voxels whose intensities are below 140. This cluster generally corresponds to the 
lowest density regions within the interior of the gas filled colon. Cluster 2 can include 

20 voxels which have intensity values in excess of 2200. These intensity values 

correspond to the enhanced stool and fluid within the colon as well as bone. Cluster 
3 can include voxels vsdth intensities in the range of about 900 to about 1 080. This 
intensity range generally represents soft tissues, such as fat and muscle, which are 
unlikely to be associated with the colon. The remaining voxels can then be grouped 

25 together as cluster 4, which are likely to be associated with the colon wall (including 
mucosa and partial volume mixtures aroimd the colon wall) as well as lung tissue and 
soft bones. 

Clusters 1 and 3 are not particularly valuable in identifying the colon 
wall and, therefore are not subject to substantial processing during image 
30 segmentation procedures for virtual colonoscopy. The voxels associated with cluster 
2 are important for segregating stool and fluid from the colon wall and are processed 
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further during the high-level feature extraction operations. Low level processing is 
concentrated on the fourth cluster, which has the highest likelihood of corresponding 
to colon tissue (step 1550). 

For each voxel in the fourth cluster, an intensity vector is generated 

5 using itself and its neighbors. The intensity vector provides an indication of the 
change in intensity in the neighborhood proximate a given voxel. The number of 
neighbor voxels which are used to establish the intensity vector is not critical, but 
involves a tradeoff between processing overhead and acciiracy. For example, a simple 
voxel intensity vector can be established with seven (7) voxels, which includes the 

0 voxel of interest, its fix>nt and back neighbors, its left and right neighbors and its top 
and bottom neighbors, all surrounding the voxel of interest on three mutually 
perpendicular axes. Figure 17 is a perspective view illustrating an exemplary intensity 
vector in the form of a 25 voxel intensity vector model, which includes the selected 
voxel 1702 as well as its first, second and third order neighbors. The selected voxel 

5 1702 is the central point of this model and is referred to as the fixed voxel. A planar 
slice of voxels, which includes 12 neighbors on the same plane as the fixed voxel, is 
referred to as the fixed slice 1704. On adjacent planes to the fixed slice are two 
nearest slices 1706, having five voxels each. Adjacent to the first nearest shces 1706 
are two second nearest slices 1708, each having a single voxel. The collection of 

0 intensity vectors for each voxel in the fourth cluster is referred to as a local vector 
series. 

Because the data set for an abdominal image generally includes more 
than 300 slice images, each with a512x512 voxel array, and each voxel having an 
associated 25 voxel local vector, it is desirable to perform feature analysis (step 1570) 

5 on the local vector series to reduce the computational burden. One such feature 
analysis is a principal component analysis (PCA), which can be applied to the local 
vector series to determine the dimension of a feature vector series and an orthogonal 
transformation matrix for the voxels of cluster 4, 

It has been found that the histogram (Figxire 16) of the CT image 

0 intensities tends to be fairly constant from patient to patient for a particular scanner, 
given equivalent preparation and scaiming parameters. Relying on this observation, 
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an orthogonal transformation matrix can be established which is a predetermined 
matrix determined by using several sets of training data acquired using the same 
scanner under similar conditions. From this data, a transformation matrix, such as a 
Karlhunen-Lo^ve (K-L) transformation, can be generated in a known manner. The 
S transformation matrix is applied to the local vector series to generate feature vector - 
series. Once in the feature- vector space domain, vector quantization techniques can 
be used to classify the feature vector series. 

An analytical, self-adaptive algorithm can be used for the classification 
of the feature vectors. In defining this algorithm, let {X^eR^a = 1,2,3,...,N} be the 

10 series of the feature vectors, where N is the number of feature vectors; K denotes the 
maximum number of classes; and T is a threshold which is adaptive to the data set. 
For each class, a representative element is generated by the algorithm. Let a,^ be a 
representative element of class k and n,^ be the number of feature vectors in that class. 
The algorithm can then be outlined as: 

15 1. Set n^=l; a^=X^;K=l; 

2. obtain the class number 'k and class parameters (aj^, n^) 
for (i=l;/<N;/-H-) 

for (i = l; j<^; J + + ) 

calculate = dist(X^,a^) ; 

20 end for 

index=arc min dj; 

m 

update class parameters: 



^index "Z 7T ^ ^ ^index* ^ index -^i ^ ' 

-"index 

^index'^index^^' 
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end if 
else 

generate new class 



end else 



end for 



3. label each feature vector to a class according to the nearest neighbor rule 
for (i = l; i<i\?; i + +) 

. for (j=l;i<K; j + + ) 
10 calculate = dist {X^ ,a^) ; 

end for 

index = arc min dj. 

label voxel / to class index. 

end for 



15 In this algorithm, dist(x,y) is the Euclidean distance between vector x 

and>' and arc min d, gives the integer j which realizes the minimum value of dy. 

The above described algorithm is dependent only on the parameters T 
and K. Hov^ever, the value of K, which relates to the number of classes within each 
voxel cluster, is not critical and can be set to a constant value, such as K'^IS. 

20 However, T, which is the vector similarity threshold, greatly influences the 

classification results. If the selected value of T is too large, only a single class will be 
generated. On the other hand, if the value of T is too small, the resulting classes will 
exhibit undesirable redundancy. By setting the value of T to be equal to the maximum 
component variance of the feature vector series, the maximum number of distinct 

25 classes results. 
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As a result of the initial classification process, each voxel within the 
selected cluster is assigned to a class (step 1570). In the exemplary case of virtual 
colonoscopy, there are several classes within cluster 4. Thus, the next task is to 
determine which of the several classes in cluster 4 corresponds to the colon wall. The 
5 first coordinate of the feature vector, which is that coordinate of the feature vector 
exhibiting the highest variance, reflects the information of the average of the 3D local 
voxel intensities. The remaining coordinates of the feature vector contain the 
information ofdirectional intensity change within the local neighbors. Because the 
colon wail voxels for the interior of the colon are generally in close proximity to the 
10 gas voxels of cliister 1, a threshold interval can be determined by data samples 
selected firom typical colon wall intensities of a typical CT data set to roughly 
distinguish colon wall voxel candidates. The particular threshold value is selected for 
each particular imaging protocol and device. This threshold interval can then applied 
to all CT data sets (acquired from the same machine, using the same imaging 
15 protocol). If the first coordinate of the representative element is located in the 

threshold interval, the corresponding class is regarded as the colon wall class and all 
voxels in that class are labeled as colon walMike voxels. 

Each colon walMike voxel is a candidate to be a colon wall voxel. 
There are three possible outcomes of not belonging to the colon wall. The first case 
20 relates to voxels which are close to the stool/liquid inside the colon. The second case 
occurs when voxels are in the lung tissue regions. The third case represents mucosa 
voxels. Clearly then, low level classification carries a degree of classification 
uncertEiinty. The causes of the low-level classification uncertainty vary. For example, 
a partial-volume effect resulting fi*om voxels containing more than one material type 
25 (i.e., fluid and colon wall) leads to the first case of uncertainty. The second and the 
third cases of uncertainty are due to both the partial volimie effect as well as the low 
contrast of CT images. To resolve the uncertainty, additional information is needed. 
Thus, a high-level feature extraction procedure is used in the present method to 
further distinguish candidates for the colon wall from other colon wall-like voxels, 
30 based on a priori anatomical knowledge of the CT images (step 1580). 
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An initial step of the high-level feature extraction procedure can be to 
eliminate the region of lung tissue from the low-level classification results. Figure 
18A is an exemplary slice image clearly illustrating the lung region 1802. The lung 
region 1 802 is identifiable as a generally contiguous three dimensional volume 
5 enclosed by colon wall-like voxels, as illustrated in Figure 18B. Given this 

characteristic, the lung region can be identified using a region growing strategy. The 
first step in this technique is to find a seed voxel within the region of growing. 
Preferably, the operator performing the CT imaging scan sets the imaging range such 
that the top most slice of the CT scan does not contain any colon voxels. As the 
1 0 interior of lung should be filled with air, the seed is provided by the low-level 

classification simply by selecting an air voxel. Once the lung region outline of Figure 
1 8B is determined, the lung volume can be removed fi:om the image slice (Figure 
18C). 

A next step in performing high-level feature extraction can be to 
IS separate the bone voxels from enhanced stool/fluid voxels in cluster 2. The bone 
tissue voxels 1 902 are generally relatively far away from the colon wall and resides 
outside the colon volume. To the contrary, the residual stool 1906 and fluid 1904 are 
enclosed inside the colon volume. Combining the a priori proximity information and 
the colon wall information obtained from the low-level classification process, a rough 
20 colon wall volume is generated. Any voxel separated by more than a predetermined 
number (e.g., 3) of voxel units from the colon wall, and outside the colon volume, 
will be labeled as bone and then removed from the image. The remaining voxels in 
cluster 2 can be assumed to represent stool and fluid within the colon volume (see 
Figures 19A-C). 

25 The voxels within the colon volume identified as stool 1 906 and fluid 

1904 can be removed from the image to generate a clean colon lumen and colon wall 
image. In general, there are two kinds of stool/fluid regions. One region type is small 
residual areas of stool 1906 attached to the colon wall. The other region type is large 
volumes of fluid 1904, which collect in basin-like colonic folds (see Figures 19A-C). 

30 The attached residual stool regions 1906 can be identified and removed 

because they are inside the rough colon volume generated during the low-level 
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classification process. The fluid 1906 in the basin-like colon fold usually has a 
horizontal surface 1908 due to the effect of gravity. Above the surface is always a gas 
region, which exhibits a very high contrast to the fluid intensity. Thus, the surface 
interface of the fluid regions can be easily marked. 
5 Using a region growing strategy, the contour of the attached stool 

regions 1906 can be outlined, and the part which is away from the colon wall volume 
can be removed. Similarly, the contour of the fluid regions 1904 can also be outlined. 
After eliminating the horizontal surfaces 1908, the colon wall contour is revealed and 
the clean colon wall is obtained. 

10 It is difficult to distinguish the mucosa voxels from the colon wall 

voxels. Even though the above three dimensional processing can remove some 
mucosa voxels, it is difficult to remove all mucosa voxels. In optical colonoscopy, 
physicians directly inspect the colonic mucosa and search for lesions based on the 
color and texture of the mucosa. In virtual colonoscopy, most mucosa voxels on the 

1 5 colon wall can be left intact in order to preserve more information. This can be very 
useful for three dimensional volume rendering. 

From the segmented colon wall volume, the inner surface, the outer 
surface and the wall itself of the colon can be extracted and viewed as a virtual object. 
This provides a distinct advantage over conventional optical colonoscopy in that the 

20 exterior wall of the colon can be examined as well as the interior wall. Furthermore, 
the colon wall and the colon lumen can be obtained separately from the segmentation. 

Because the colon is substantially evacuated prior to imaging, a 
commonly encountered problem is that the colon lumen collapses in spots. While the 
inflation of the colon with compressed gas, such as air or CO2, reduces the frequency 

25 of collapsed regions, such areas still occur. In performing a virtual colonoscopy, it is 
desirable to automatically maintain a flight path through the collapsed regions and it is 
also desirable to use the scanned image data to at least partially recreate the colon 
lumen in the collapsed regions. Since the above described image segmentation 
methods effectively derive both the interior and exterior of the colon wall, this 

30 information can be used to enhance the generation of the fly path through the 
collapsed regions. 
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In extending the flight path through collapsed regions of the colon or 
expanding a collapsed region of the colon, the first step is to detect a collapsed region. 
Using the premise that the grayscale values of the image data from around the outside 
of the colon wall change much more dramatically than the greyscale values within the 
5 colon wall itself, as well as in other regions such as fat, muscle and other kinds of 
tissue, an entropy analysis can be used to detect areas of colon collapse. 

The degree of change in greyscale value, for example along the 
centerline, can be expressed and measured by an entropy value. To calculate an 
entropy value, voxels on the outer surface of the colon wall are selected. Such points 
10 are identified fi'om the above described image segmentation techniques. A 5x5x5 
cubic window can be applied to the pixels, centered on the pixel of interest. Prior to 
calculating the entropy value, a smaller (3x3x3) window can be applied to the pixels 
of interest in order to filter out noise fi^om the image data. The entropy value of a 
selected window about the pixel can then be determined by the equation: 

E=Y^C{i)ln{C{i)) 



1 5 where £ is the entropy and C(i) is the number of points in the window with the 

grayscale of i (i=0,l,2,. . ., 255). The calculated entropy values for each window are 
then compared against a predetemiined threshold value. For regions of air, the 
entropy values will be fairly low, when compared to regions of tissue. Therefore, 
along the centerline of the colon lumen, when the entropy values increase and exceed 

20 the predetermined threshold value, a collapsed region is indicated. The exact value of 
the threshold is not critical and will depend in part on the imaging protocol and 
particulars of the imaging device. 

Once a collapsed region is detected, the previously determined 
centerline flight path can be extended through the region by piercing through the 

25 center of the collapse with a one voxel wide navigation line. 

In addition to automatically continuing the flight path of the virtual 
camera through the colon lumen, the region of colon collapse can be virtually opened 
using a physical modeling technique to recover some of the properties of the collapsed 
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region. In this technique, a model of the physical properties of the colon wall is 
developed. From this model, parameters of motion, mass density, damping density, 
stretching and bending coefficients are estimated for a Lagrange equation. Then, an 
expanding force model (i.e., gas or fluid, such as air, pumped into the colon) is 

5 formulated and applied in accordance with the elastic properties of the colon, as 

defined by the Lagrange equation, such that the collapsed region of the colon image is 
restored to its natural shape. 

To model the colon, a finite-element model can be applied to the 
collapsed or obstructed regions of the colon lumen. This can be performed by 

0 sampling the elements in a regular grid, such as an 8 voxel brick, and then applying 
traditional volume rendering techniques. Alternatively, an irregular volume 
representation approach, such as tetrahedrons can be applied to the collapsed regions. 

In applying the external force (air pumping) model to the colon model, 
the magnitude of the external force is first determined to properly separate the 

S collapsed colon wall regions. A three dimensional growing model can be used to 
trace the internal and external colon wall surfaces in a parallel manner. The 
respective surfaces are marked from a starting point at the collapsed region to a 
growing source point, and the force model is applied to expand the surfaces in a like 
and natural manner. The region between the internal and external surfaces, i.e., the 

0 colon wall, are classified as sharing regions. The external repulsive force model is 
applied to these sharing regions to separate and expand the collapsed colon wall 
segments in a natural manner. 

To more clearly visualize the features of a virtual object, such as the 
colon, which is subjected to virtual examination, it is advantageous to provide a 

5 rendering of the various textures of the object. Such textures, which can be observed 
in the color images presented during optical colonoscopy, are often lost in the black 
and white, grey scale images provided by the CT image data. Thus a system and 
method for texture imaging during virtual examination is required. 

Figure 20 is a flow chart depicting a present method for generating 

0 virtual objects having a texture component. The purpose of this method is to map 
textures obtained by optical colonoscopy images in the red-green-blue (RGB) color 
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space, as for example from the Visible Human, onto the gray scale monochrome CT 
image data used to generate virtual objects. The optical colonoscopy images are 
acquired by conventional digital image acquisition techniques, such as by a digital 
"frame grabber" 1429 which receives analog optical images from a camera, such as a 
5 video camera, and converts the image to digital data which can be provided to CPU 
1423 via interface port 143 1 (Figure 14). The first step in this process is to segment 
the CT image data (step 2010). The above described image segmentation techniques 
can be applied to choose intensity thresholds in the grey scale image to classify the CT 
image data into various tissue types, such as bone, colon wall, tissue, air, and the like. 

10 In addition to performing image segmentation on the CT image data, 

the texture features of the optical image need to be extracted from the optical image 
data (step 2020). To do this, a gausian filter can be apphed to the optical image data 
followed by sub-sampling to decompose the data into a multiresolutional pyramid. A 
laplacian filter and steerable filter can also be applied to the multiresolutional pyramid 

1 5 to obtain oriented and non-oriented features of the data. While this method is 

effective at extracting and capturing the texture features, the implementation of this 
approach requires a large amount of memory and processing power. 

An alternative approach to extracting the texture features from the 
optical image is to utilize a wavelet transform. However, while wavelet 

20 transformations are generally computationally efficient, conventional wavelet 

transforms are limited in that they only capture features with orientations parallel to 
the axes and cannot be appUed directly to a region of interest. To overcome these 
hmitations, a non-separable filter can be employed. For example, a lifting scheme can 
be employed to build filter banks for wavelets transform in any dimension using a two 

25 step, prediction and updating approach. Such filter banks can be synthesized by the 
Boor-Rom algorithm for multidimensional polynomial interpolation. 

After the textural features are extracted from the optical image data, 
models must be generated to describe these features (step 2030). This can be 
performed, for example, by using a non-parametric multi-scale statistical model which 

30 is based on estimating and manipulating the entropy of non-Gausian distributions 
attributable to the natural textures. 
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Once texture models are generated fiom the optical image data, texture 
matching must be performed to correlate these models to the segmented CT image 
data (step 20S0). In regions of the CT image data where the texture is continuous, 
corresponding classes of texture are easily matched. However, in boundary regions 
5 between two or more texture regions, the process is more complex. Segmentation of 
the CT data aroimd a boundary region often leads to data which is fuzzy, i.e., the 
results reflect a percentage of texture from each material or tissue and vary depending 
on the various weighting of each. The weighting percentage can be used to set the 
. importance of matching criteria. 

10 In the case of the non-parametric multi-scale statistical model, the 

cross entropy or a Kullback-Leiber divergence algorithm can be used to measure the 
distribution of different textures in a boimdary region. 

After texture matching, texture synthesis is performed on the CT image 
data (step 2050). This is done by fusing the textures from the optical image data in to 

15 the CT image data. For isotropic texture patterns, such as presented by bone, the 

texture can be sampled directly from the optical data to the segmented CT image data. 
For anisotropic texture regions, such as colon mucosa, a muitiresolution sampling 
procedure is preferred. In this process, selective re-sampling for homogenous and 
heterogenous regions is employed. 

20 Alternatively, pseudocolor texture can be created directly from the CT 

data. For each voxel, multiple CT values, comprising a local area neighborhood, can 
be evaluated to determine a pseudocolor for the given voxel. For each voxel the local 
neighborhood consists of the voxels that are within some given distance of the center 
voxel. For example a 5x5x5 voxel cubic shaped region, or the double pyramid which 

25 represents all voxels within 3 units measured by Manhattan distance. This vector of 
scalar values is then evaluated to map to a color to be displayed for this voxel during 
subsequent volume rendering. The evaluation of the local neighborhood vector of 
values can compute such things as local curvature, homo/heterogeneity, or other 
geometric or spatial functions. 



30 



Volume Rendering 
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In addition to image segmentation and texture mapping described 
above, volume rendering techniques can be used in connection with virtual 
colonoscopy procedures to fiirther enhance the fidelity of the resulting image. Figure 
21 illustrates a perspective volume ray-casting method which can be used for volume 
5 rendering in accordance with the present invention. From a selected virtual 
viewpoint, e.g., camera position, such as within the colon lumen, rays are cast 
through each of the proximate image pixels (step 2100). For each ray, the first 
sampling point is set as the current im^e pixel along the ray (step 2110). The 
distance (d) between the current sampling point and the nearest colon wall is then 
0 determined (step 2120). The current distance (d) is compared to a predetermined 
sampling interval (i) (step 2130). If the distance (d) is greater than the sampling 
interval (i) then no sampling occurs and the next sampling point along the ray is 
determined by jumping the distance d along the ray (step 2140). If the distance is less 
than or equal to the sampling interval (i) then conventional sampling is performed on 
this point (step 2150) and the next sampling point is selected in accordance with the 
sampling interval (i) (step 2160). For example, trilinear interpolation between the 
density values of 8 neighboring voxels can be performed to determine the new density 
value at the sampling point. 

The method of Figure 21 effectively accelerates ray-casting because a 
space leaping technique is used to quickly skip over empty space along the ray of the 
image plane to the colon wall.. In this method, a distance fi'om a sample point to the 
nearest colon wall is determined along each ray. If the distance is larger than a 
predetermined sampling interval (i), a jump to the next sampling point along the ray is 
performed. Since the closest distance information is already available fi-om the 
potential field which is used for virtual camera control, no additional distance coding 
calculations are required. In this case, neither surface rendering nor Z-buffer 
transform is required, which results in savings in preprocessing time and memory 
spaced 

Alternatively, a space leaping method can derive distance information 
for each ray firom the Z-buffer of the corresponding surface rendering image. If the 
surface rendering image and volume rendering image will both be generated, this 
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approach provides minimal processing overhead burden as the Z-buffer information is 
provided as a result of the surface rendering methods. Thus, this form of space leaping 
method only requires additional processing to perform a depth transformation from 

♦ 

the image space domain to the w^orld space domain. 
5 For those regions along the ray where the distance (d) was traversed in 

step 2140, the region along the ray corresponds to open space and can be assigned a 
value according to an open space transfer function. Typically, open space will have 
no contribution on the final pixel value. For each point where sampling takes place, 
one or more defined transfer functions can be assigned to map different ranges of 
10 sample values of the original volume data to different colors and opacities and 
possibly other displayable parameters. For example, four independent transfer 
functions have been used to determine different material by mapping ranges of CT 
density values into specified colors of red, green, blue and opacity, each in the range 
of0to255, 

15 Virtual Biopsy 

The above described techniques can also form the basis of a system for 
performing virtual electronic biopsy of a region being examined to effect a flexible 
and non-invasive biopsy. As noted above, volume rendering techniques use one or 
more defined transfer functions to map different ranges of sample values of the 

20 original volume data to different colors, opacities and other displayable parameters for 
navigation and viewing. During navigation, the selected transfer function generally 
assigns maximimi opacity to the colon wall such that the outer surface is easily ' 
viewed. Once a suspicious area is detected during virtual examination, the physician 
can interactively change the transfer function assigned during the volume rendering 

25 procedure such that the outer surface being viewed becomes substantially transparent, 
allowing the region information to be composited and thus the interior structure of the 
region to be viewed. Using a number of predetermined transfer functions, the 
suspicious area can be viewed at a number of different depths, with varying degrees of 
opacity assigned throughout the process. 
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Polyp Detection 

The present system and methods can be used to perform automated 
polyp detection. With reference to Figure 13, polyps 1305, which occur, for example, 
within the colon, generally take the form of small convex hill-like structures 
5 extending from the colon wall 1301. This geometry is distinct from the fold of the 
colon wall. Thus, a differential geometry model can be used to detect such polyps on 
the colon wall. 

- The surface of the colon lumen can be represented as a continuously 
second differentiable surface in three dimensional Euclidean space, such as by using a 

1 0 C-2 smoothness surface model. Such a model is described in "Modem Geometry 
Methods and Applications" by B.A, Dubrovin et al, published by Springer- Verlag 
1994, which is hereby incorporated by reference in its entirety. In this model, each 
voxel on the surface of the colon has an associated geometrical feature which has a 
Gauss curvature, referred to as Gauss curvature fields. A convex hill on the surface, 

IS which may be indicative of a polyp, possesses a unique local featiure in the Gauss 
curvature fields. Accordingly, by searching the Gauss curvature fields for specific 
local features, polyps can be detected. Once detected, the suspected polyps can be 
highlighted and thus brought to the attention of the physician where the physician can 
measure the suspected polyp and use the above described virtual biopsy methods to 

20 further investigate the suspicious region. 



Central Flv-Path Generation 

In the case of virtual colonoscopy, determining a proper navigation 
line, or fly-path, through the colon lumen is an hnportant aspect of the described 
systems and methods. While certain techniques for determining the fly-path of the 

25 virtual camera model were discussed with respect to Figures 4-8, Figure 22 illustrates 
an alternate method of generating the central fly-path through the colon lumen. After 
the colon wall is identified, such as by the image segmentation methods described 
herein, a volume shrinking algorithm can be employed to emphasize the trend of the 
colon lumen and reduce subsequent searching time within the lumen volume (step 

30 2310). 
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Figure 23 further illustrates the steps of an exemplary volume 
shrinking algorithm, which is based on a multiresolution analysis model. In this 
procedure, the three dimensional volume is represented by a stack of binary images 
which have the same matrix size (step 23 10). Collectively, these images form a 
5 binary data set. A discrete wavelet transformation can be applied to the binary data 
set which results in a number of sub-data sets representing different time-frequency 
components of the binary data set (step 2320). For example, the discrete wavelet 
transformation may yield eight (8) sub-data sets. The sub-data sets are compared 
against predetermined threshold values such that the lowest frequency component is 

10 identified (step 2330). This component forms the binary data set for subsequent 

discrete wavelet transformation and thresholding steps, which are recursively applied 
in a multi-resolution structure (step 2340). In the case of virtual colonoscopy, the 
discrete wavelet transformation and associated thresholding can be applied three times 
recursively on the subsequent sub-dataset that represents the lowest frequency 

1 S component (a 3-level multi-resolution decomposition). 

Returning to Figure 22, from the reduced colon volume model, a 
distance map technique can be employed to generate a minimum distance path 
between the two ends of the colon, e.g., from the rectum to the cecum (step 2215). 
The resulting path preserves the global trend information of the colon lumen, but 

20 ignores the trends exhibited by local folds. Control points within the global colon can 
then be determined by mapping the minimum distance path back to the original data 
space (Step 2220). For example, in the case of a 3-level multi-resolution 
decomposition, the reduced volume is three times smaller than the original volume 
and an affine transformation, which is well known, can be used to map the reduced 

25 volume model exactly back to the original scale volume. The minimum distance path 
of the reduced value can also be mapped back into the original scale volume as a 
series of points, which can be used as the control points within the colon. 

The preferred fly path is one which is on the centerline of the colon 
lumen. However, the initial control points may not be exactly located in the center of 

30 the colon lumen. Thus, the initial control points can be centered, such as by the use of 
a bi-section plane algorithm (step 2230). For example, at each selected control point. 
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a bi-section plane can be defined as a plane normal to the trend direction and cutting 
across the colon lumen. A centralization algorithm, such as a maximum disk 
algorithm, can then be performed on each bi-section plane. Such an algorithm is 
discussed in the article "On the Generation of Skeletons from Discrete Euclidean 

* 

5 . Distance Maps" by Ge et al., ffiEE Transactions on PAMI, Vol. 1 8, pp. 1055-1 066, 
1 996 which is hereby incorporated by reference. 

Once the control points are centralized, the flight path can be 
determined by interpolating a line connecting these points (step 2240). hi the case of 
virtual colonoscopy, it is desirable that the interpolated flight path take the form of a 

0 smooth curve which is substantially centered within the colon lumen. A constrained 
cubic B-spline interpolation algorithm based on Serret-Frenet Theorem in differential 
geometry theory can be used to establish a suitable smooth curved flight path, such as 
is described in "Numerical Recipes in C: The Art of Scientific Computing," by Press 
et al,, Second Edition, Cambridge University Press, 1992. 

5 The pictorial representation of a segmented colon lumen in Figures 24 

and the flow chart of Figure 25 s^et forth yet another alternate fly-path generation 
method in accordance with the present invention. In this alternate method, the 
representation of the colon lumen 2400 is first partitioned into a number of segments 
2402 a-g along the length of the lumen 2400 (step 2500). From within each segment 

0 2402 a representative point is selected 2404 a-g (step 2520). Each representative 
point 2404 a-g is then centered with respect to the colon wall (step 2530), such as by 
the use of a physically-based deformable model which is used to push the points to the 
center of the respective segment. After the representative points are centered, the 
points are sequentially joined to establish the center-line fly-path for the virtual 

15 camera model (step 2540). If the segments are sufficiently small in length, the 
centered points can be connected with straight line segments 2406 a-f However, . 
when linear curve fitting techniques are applied to join the centered points, a 
smoother, continuous flight path is established. 

Each of the foregoing methods can be implemented using a system as 

10 illustrated in Figure 14, with appropriate software being provided to control the 
operation of CPU 1409 and CPU 1423. 
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An alternate hardware embodiment, suitable for deployment on a 
personal computer, is illustrated in Figure 26. The system includes a processor 2600 
which should take the form of a high speed, multitasking processor such as a Pentium 
in processor operating at a clock speed in excess of 400 MHZ. The processor 2600 is 

5 coupled to a conventional bus structure 2620 which provides for high speed parallel 
/ data transfer. Also coupled to the bus structure 2620 are main memory 2630, a 
graphics board 2640, and a volume rendering board 2650. The graphics board 2640 is 
preferably one which can perform texture mapping, such as the Diamond Viper v770 
Ultra manufactured by Diamond Multimedia Systems. The volume rendering board 

) 2650 can take the form of the VoluraePro board from Mitsubishi Electric, which is 
based on U.S. Patent Nos. 5,760,781 and 5,847,71 1, which are hereby incorporated by 
reference. A display device 2645, such as a conventional SVGA or RGB monitor, is 
operatively coupled to the graphics board 2640 for displaying the image data. A 
scanner interface board 2660 is also provided for receiving data from an imaging 

> scanner, such as an MRI or CT scanner, and transmitting such data to the bus structure 
2620. The scanner interface board 2660 may be an application specific interface 
product for a selected imaging scanner or can take the form of a general purpose 
input/output card. The PC based system 2600 will generally include an I/O interface 
2670 for coupling I/O devices 2680, such as a keyboard, digital pointer (e.g., mouse) 

) and the like to the processor 2620. Alternatively, the I/O interface can be coupled to 
the processor 2620 via the bus 2620. 

In the case of three dimensional imaging, including texture synthesis 
and volume rendering, numerous data handling and processing operations are 
required. For large datasets, such as those represented by the colon lumen and its 

5 surrounding area, such processing can be very time consuming and memory intense. 
However, using the topology of Figure 26 in accordance with the processing method 
illustrated in the flow chart of Figure 27, such operations can be performed on a 
relatively low cost personal computer (PC). Imaging data is received by the processor 
2620 and stored in main memory 2630 via the scanner interface board 2660 and bus 

Q structure 2620. This image data (pixels) is converted into a volxune element (voxel) 
representation (step 271 0). The volume representation, which is stored in main 
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memory 2630, is partitioned, for example into slices, such as along a major volume 
axis or other portions of the region being imaged (step 2720). The volume partitions 
are then transferred to the volume rendering board and temporarily stored in volume 
rendering memory 2655 for volume rendering operations (step 2730). The use of 

5 locally resident volume rendering memory 2655 provides for enhanced speed in 
volume rendering as data need not be exchanged over the bus 2620 during rendering 
of each slice of the total volume. Once volume rendering is complete for the slice, the 
rendered data is transferred back to main memory 2630 or the graphics board 2640 in 
a sequential buffer (step 2740). After all slices of interest have been subjected to 

10 rendering, the contents of the sequential buffer are processed by the graphics board 
2640 for display on tiie display unit 2645 (step 2750). 

Multi-scan Based Virtual Examination 

The techniques discussed above generally perform virtual imaging 
based on a dataset acquired from a single magnetic resonance imaging (MRI) or 

1 5 computed tomography (CT) scan. However, the techniques discussed above are also 
useful for performing virtual examination of a region using multiple scans of a region. 
By using multiple scans of a region, improved imaging of regions of pathology can be 
achieved and motion artifacts can be reduced. One such application of interest is in 
performing virtual cystoscopy to screen a patient for possible polyps or cancer of the 

20 bladder. 

Figure 28 is a flow chart which illustrates a method of employing 
multiple MRI scans to perform virtual examination of an object, such as virtual 
cystoscopy. Unlike CT images, where the bladder wall can be difficult to distinguish 
from urine, in MRI images, urine can be used as a natural contrast agent to delmeate 

4 

25 the inner bladder wall. To this end, a pre-image scan protocol is employed (step 2805). 
Approximately 14 hour prior to the first of four MRI scans, the patient is requested to 
empty the bladder and then consume one cup of water. After approximately hour, 
the patient is subjected to the first of four MRI scans of the bladder region (step 28 1 0). 
The first scan, with the bladder full and distended, follows protocol for Tl -weighted 

30 transverse imaging. For example, when using the Picker scanner referenced above, a 
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KJELL FASTER protocol using a 256x256 matrix size, a 38 cm field of view (FOV), a 
1.5 mm slice thickness (no gap), a 3 ms TE, a 9 ms TR, a 30 degree flip angle and one 
scan average can be used. Of course, these parameters tend to be scanner specific and 
various changes in the parameters can be used with acceptable results. 
5 With the bladder still full, the patient is subjected to a second MRI scan, 

scan 2 (step 2815). The second scan follows a protocol for Tl-Aveighted coronal 
imaging, such as the Picker KJELL FASTER protocol with a 256x256 matrix size, a 
38 cm field of view (FOV), a 1 .5 mm slice thickness (no gap), a 3 ms TE, a 9 ms TR, a 
30 degree flip angle and a two-scan average. 

1 0 The two image scans described above are teiken along orthongal axes 

with respect to one another. The advantage of this is that regions of significant motion 
artifacts in one scan, generally correspond to regions of minimal motion artifacts in the 
orthogonal scan. Accordingly, by taking a first scan in the transverse direction and a 
second scan in the coronal direction, the image scans can be registered and motion 

1 5 artifacts in the data set can be identified and compensated for. 

After the scan 2, the patient is asked to relieve the bladder and is then 
subjected to two additional MRI scans. The third scan (step 2820) follows the same 
imaging protocol as the first scan (transverse imaging). The fourth scan (step 2825) 
follows the same imaging protocol as the second scan (coronal imaging). 

20 The image scans can be acquired using a Picker 1 .5 T Edge whole-body 

a 

scanner. Although a T2 imaging protocol can be used, a Tl imaging protocol is 
preferred for virtual cystoscopy because this protocol provides improved delineation 
between fat and mine and offers a shorter acquisition period. Alternatively, the image 
scans can take the form of computed tomography or ultrasoimd imaging scans using 

25 suitable contrast agents and protocols for these imaging techniques. 

During the first two scans (scan 1 and scan 2), the bladder is distended 
and the bladder wall is relatively thin. In this case, physiologically altered locations, 
such as tumors, may thin at a different rate as compared to the unaltered bladder wall 
and may become more apparent imder these conditions. During the third and fourth 

30 scans, the bladder is substantially empty and the bladder wall is thicker. With a thicker 
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wall, a more pronounced image contrast may result between normal tissue of the 
bladder wall and that of physiologically altered tissue. 

After the four scans are acquired, the four corresponding datasets can 
then be individually processed. Initially, each scan data set is preferably subjected to 

S image segmentation, as discussed above (step 2830), such as in connection with Figure 
15. During image segmentation, the voxels of the four datasets are classified into a 
number of categories, such as bladder wall, urine, fat, boundary, etc. The classification 
is based on the local intensity vectors of the voxels. Once the voxels are classified, the 
interior of the bladder lumen can be identified using a region growing algorithm 

0 begiiuiing with a seed voxel selected firom within the bladder volume, such as by 
selecting an air voxel or urine voxel. 

Prior to clinical analysis of the segmented volume data sets, registration 
of the four data sets to a common coordinate system is performed (step 2835). Because 
the shape of the bladder varies from scan to scan, an exact voxel-voxel registration is 

5 not of practical value. Instead a flexible registration process is preferred. In the 

present flexible registration process, for each volume of interest (volume rendered for 
each corresponding scan) the center of the volume is determined, such as by averaging 
the three coordinates of all the voxels in the volume. 

A Cartesian coordinate system can then be constructed with the origin 

0 of the system located at the center point of the volume. The axes of the system can 
then be oriented in a number of ways. A suitable selection of orientation corresponds 
to the orientation of the natural human body, e.g., with the Z-Axis running along the 
height of the body (e.g., from toe to head) the Y-axis oriented firom back to fi'ont and 
the X-axis running laterally ( e.g., firom left to right). The units of length in this 

5 coordinate system can be conveniently set to an arbitrary unit of one voxel length, the 
absolute magnitude of which will vary based on acquisition properties for the MRI 
scans. So long as the same pixel spacing is used in all scans to acquire all four data 
sets, this will result in a uniform value for each of the four data sets. 

After registration, the images from the four data sets can be viewed 

0 individually or simultaneously (step 2845). An exemplary display window is 

illustrated in Figures 29 and 30. Referring to Figure 29, the display is partitioned into 
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four sub windows 2905, 2910, 2915, 2920 which correspond to scan 1, scan 2, scan 3 
and scan 4, respectively. A control panel section 2925 can also be provided on a 
portion of the display to establish a graphical user interface (GUI) to offer display and 
navigation functions to the user. As an operator navigates in one of the image sub 

5 windows, such as magnifying the view, the corresponding operation preferably takes 
place in the other sub window views as well. A user can also select one of the views 
for expansion to a single window display. 

To reduce the amount of data which is simultaneously processed, the 
data sets can be partitioned, such as into 8 parts or octants (step 2840). This can be 

0 performed in a number of ways. For example, with reference to the Cartesian 
coordinate system illustrated in Figure 29, the data can be portioned into the eight 
regions of the coordinate system: (1) X, Y, Z; (2) X, -Y, Z; (3) X, Y, -Z; (4) X, -Y, -Z; 
(5) -X, Y, Z; (6) -X, - Y, Z; (7) -X, Y, -Z; and (8) -X, -Y, -Z. 

Figure 29 illustrates four views of the outside of the bladder lumen 

5 taken from each of the four scans. Figure 30 illustrates fours views of a portion of the 
interior of the bladder lumen also taken from each of the four scans. 

Multi-Resolution Imaging and Virtual Larvngoscopv 

The systems and methods described herein can be adapted and applied 

0 to perform multiresolution unaging which is well suited for virtual laryngoscopy. 
Figure 3 1 is a flow chart illustrating a method for performing virtuallaryngoscopy. 
First, an imaging scan of the region of a patient's larynx is acquired (step 3105). This 
can be performed using computed tomography (CT) or magnetic resonance imaging 
(MRI) techniques. However, because the CT scan in this region offers significantly 

5 faster acquisition time (30 seconds versus over 7 minutes for MRI) and higher 

resolution (0.3nuii cubic voxel compared to 1mm cubic voxel for MRI), the CT scan is 
preferred. To acquire the CT scan data a GE/CTI spiral scan CT scanner can be used. 
A suitable scan protocol is 120 keV, 200 ma, 512x512 matrix size, 15 cm field-of-view 
and 3mm/2.0:l pitch. The scan is completed in approximately 30 seconds and results 

0 in 35 1 image slices of 0.3mm thickness and results in 0.3mm cubic voxels. 
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Image segmentation can be used to classify voxels into a number of 
categories (step 3110). In this operation, a modified self-adaptive on-line vector 
quantization (SOVQ) algorithm can be used. In such a case, the algorithm analyzes 
each voxel with respect to neighbors of up to the third order to determine local density 
5 features. Each voxel in the acquired dataset has an associated local density vector. By 
transforming the local density vectors using the Karhunen-Loeve (K-L) transform, 
feature vectors for the voxels in the volume image can be obtained. Based upon the 
feature vectors, the voxels can be classified and labeled. Voxel classification is 
dependent in part on the choice of a local voxel density vector and one preset 
parameter, referred to as the maximum cluster number (MCN). The MCN sets the 
number of voxel classifications that will be applied to the dataset. In the case of the 
CT images, the human eye can discern four (4) distinguishable tissue/material types. 
An MCN value of S is suitable in this case. For an MRI image, the human eye can 
differentiate among 6 different tissue types, and an MCN value of 7 can be used. 

As part of the image segmentation process, an expanded data set is 
generated by interpolation between the measured data points. For example, prior to 
employing the SOVQ algorithm, a first order Lagrange interpolation can be applied to 
each slice in the dataset. This expands the 256x256 matrix size of the original slices of 
the data set to a 512 x 512 matrix size. In addition, inter-slice interpolation can be 
performed to fiirther expand the dataset between actual slices. The interpolated dataset 
is referred to as the enlarged dataset. In addition to generating an enlarged dataset, the 
interpolation process also suppresses noise and reduces the partial-volume effect, as 
the interpolation process has a low-pass filtering effect on the data. 

Using a two dimensional viewing tool, a seed voxel can be selected 
within the larynx lumen and a growing algorithm applied to extract the larynx volume 
firom the dataset (step 3115). In those regions of the larynx where there may be several 
unconnected volume regions, multiple seed points can be selected. 

With the larynx volume identified and the voxels of the regions 
classified through image segmentation, the next task is to manage the data in a manner 
which allows efficient navigation and viewing of the virtual larynx. In this case, a 
level-of-detail (LOD) approach is adopted and modified for use in the present method. 
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In this LOD method, a reduced dataset is generated from the enlarged data set. For 
example, the 512x512x256 enlarged dataset can be reduced to a 64 x 64 x 32 reduced 
volume dataset using a multi-resolution decomposition with three levels of 
thresholding (step 3120). Next, polygons used to render the volume images in both the 
5 enlarged and reduced volume datasets can be extracted. A traditional Marching Cubes 
method can be used to extract polygons to fit the surface of the larynx. 

One problem encountered in the prior art is managing the large number 
of polygons required to generate the three dimensional image for the enlarged dataset. 
This problem is solved in the present method by organizing the enlarged dataset in a 

1 0 Binary Space Partitioning (BSP) tree data structure (step 3130). The original image 
volume is selected as the root of the tree. The space is then partitioned into two 
subspaces containing an approximately equal number of polygons. This subdivision 
process is iteratively repeated until the number of polygons in each resulting subspace 
is below a threshold value. The threshold value can vary based on system performance 

1 5 and application requirements. The last resulting subspaces are referred to as leaf nodes 
of the tree. Once the subdivision process is complete, all of the voxels of the expanded 
dataset are stored in the leaf nodes of the BSP tree. 

During navigation or viewing, polygon culling can be applied by first 
removing those leaf nodes that are completely outside the field-of-view from current 

20 processing operations. The remaining polygons are recalled from the BSP tree, 
ordered and rendered in those spaces which were not culled. Thus, the BSP tree 
provides an effective tool for selecting a relevant portion of the dataset for a particular 
navigation or display operation. 

The enlarged and reduced datasets are cooperatively used in a two level 

25 LOD rendering mode. If a user is interacting with the object (step 3135), such as 

rotating, shifting or effecting other changes in the field of view, the polygons from the 
reduced dataset (64-sized) are rendered (step 3 140). Because of the significantly lower 
number of polygons involved, interaction with the reduced dataset volume can be 
performed faster and with less processing overhead. The tradeoff for the increased 

§ 

30 speed is reduced image resolution. If there is no interaction from the user after a 
predetermined time period, the polygons of the enlarged dataset (5 1 2-sized) are 
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selected from the BSP tree and are rendered to provide a high resolution image of the 
current field of view (step 3 145). 

Virtual Angiography 

The techniques for virtual imaging and navigation can also be adapted 
5 and applied to virtual angiography. This technique can be used for detection and 
measurement of various abnormalities and disease of the circulatory system. 

One such application of virtual angiography is the detection of 
abdominal aortic aneurysms, which generally start as small enlargements of the aortic 
vessel and exhibit a greater risk to rupture with increasing size of the aneurysm. 

1 0 Previously, the only effective method of treatment was open surgery, placing a graft 
within the aorta at the level of the aneurysm. However, this procedure has a high 
degree of associated morbidity and mortality. Recently developed per cutaneous 
placed aortic stent graft techniques have a significantly lower complication rate. 
Virtual angiography is an effective method to help plan these less invasive procedures 

1 5 and can also be an effective tool for detecting the presence of an aneurysm and tracking 
the growth of an aneurysm to determine if and when surgery is indicated. 

Figure 32 in a flow chart which provides an overview of the present 
virtual angiography method. In performing a virtual angiography, an image scan of the 
vessel, such as the aorta must be acquired (step 3205). Various imaging techniques 

20 can be used, such as Computed Tomography (CT), Magnetic Resonance Imaging 

(MRI) and ultrasound. However, an aortic CT scan is generally preferred because of 
the contrast between blood, soft tissue and calcium deposits which results in the CT 
image. 

Once an image scan data set is acquired, image segmentation techniques 
25 are then applied to the data set to classify the voxels of the dataset into a number of 
categories (step 3210). The image segmentation techniques described above, such as 
in cormection with Figure 15, are generally applicable. In this case, the various feature 
vector values of the voxels will be grouped according to categories such as blood, soft 
tissue and calcium deposits. Using a blood voxel as a seed, a region growing algorithm 
30 can be used to determine the volume and extent of the aortic lumen. 
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In the CT image, an aneurysm has image features which closely 
resemble the neighboring soft tissue. As a result, the full contour of the aneurysm can 
be difficult to establish. However, regions with calcium deposits offer significant 
contrast on the CT scan and can be used to identify portions of the aneurysm, such as 
5 the endpoints of the aneurysm on the vessel wall (step 3215). 

After a portion of an aneurysm is detected, one or more closing surfaces 
can be generated to define an estimation of the aneurysm's contour (step 3220). A 
convex closing surface can be established using a non-uniform, non-rational B-spline 
to generate a surface which runs through or near the points of the aneurysm which were 
10 identified. 

After the closing surface is generated, the volvune of the aneurysm can 
be estimated (step 3225). One method for estimatmg the volume is to count the 
number of voxels which are enclosed by the estimated closing surface. In addition, 
within the volume of the aneurysm, the centerline along the direction of blood flow can 

15 be determined by using a distance transform technique. Continuous local coordinate 
systems can then be established along this centerline and the diameter of the aneurysm 
determined. Virtual navigation can take place along this centerline, in a manner 
consistent with that described above for navigating through a lumen, such as the colon. 

Referring to Figures 33A-C, the described method of virtual 

20 angiography can be used to assist in the generation and placement of a stent graft to 
bypass an abdominal aortic aneurysm. Figure 33A illustrates a simplified diagram of 
an abdominal aortic aneurysm located below the renal arteries and above the 
bifurcation of the aorta. Because of variations from patient to patient in the specific 
anatomy of the aorta and the size and location of an abdominal aortic aneurysm therein, 

25 when a stent graft is to be used to bypass an aneurysm, the graft must be designed and 
built to specifically fit the particular aortic segment. As illustrated in Figure 333, this 
can require identifying the length of the required graft, the diameter at the points of 
interface on each end of the bypassed region, the angles of interface, among other 
variables. If the aneurysm is located near an arterial branch, the size and angles of the 

30 bifurcated ends of a bifurcated stent graft must also be determined, as illustrated in 
Figures 33B and 33C. 
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To date, such measurements have been performed through invasive 
calibrated angiograms using a catheter inserted into the aorta from an insertion made at 
the level of the groin region, rapid injection of a large amount of iodinated contrast and 
rapid radiographic imaging. This technique can be supplemented and perhaps 
5 supplanted using the present virtual angiography techniques, which can resolve such 
distances and angles using virtual navigation using centerlines constructed through the 
branches of the aortic lumen. In addition, the virtual angioscopy can be used to 
perform a virtual biopsy of the region where a stent graft may be inserted. This allows 
the operator to view beneath the arterial surface and examine the region for thrombus 
10 deposits, calcification or other factors which would contra-indicate the use of a stent 
graft procedure. 

Another application of virtual angiography is the imaging, examination 
and navigation through the carotid arteries which supply blood flow to the brain. The 
techniques described herein with respect to virtual endoscopy are fiilly applicable in 

1 5 the case of blood vessels. For example, the vessels of interest are extracted from the 
acquired image data using image-segmentation techniques. Next, a navigation flight 
path can be established through die vessel(s). Preferably, potential fields are built up 
within the vessel for use in navigation. As with other organs, such as the colon, a 
volume-rendered model of the vessels of interest can be generated. Using the flight 

20 path and potential fields to navigate through the interior of the volimie rendered blood 
vessel lumen, abnormalities such as vessel narrowing and plaque build up can be 
observed. Li addition, the techniques discussed regarding virtual biopsy can be applied 
in this context to evaluate vessel wall and characterize build up on the wall surface, 
such as plaque. 



25 

Tree Branch Searching for Virtual Endoscopy 

Path planning for virtual navigation through a hollow organ or object is 
an important task. Various techniques have been discussed, such as fly-path 
generation, to achieve this goal. As the geometry of the object being studied becomes 
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more complex, such as presenting a multi-branch structure, the task of path planning 
becomes even more complex. It is desirable to determine not only the center line of a 
primary lumen, but also to identify and locate any branches extending from the primary 
lumen. Common examples of organs having a complex branch structure include the 
S main airway and lungs, the cardiovascular system and, because of the presence of 
haustral folds, the colon. Each organ or object generally presents specific challenges 
for defming a path, or skeleton, for the object. However, a generalized technique for 
generating such a skeleton is illustrated in the flow chart of Figure 34. 

Referring to Figure 34, an imaging scan of the region of interest, such as 

1 0 a computed tomography (CT) or Magnetic Resonance Imaging (MRI) scan, is acquired 
(step 3405). As discussed above, the imaging scan is transformed into a three 
dimensional volume of the region by stacking the binary images of the imaging scan 
and defining three dimensional volume units, or voxels, from these stacked images 
(step 3410). Depending on the volume and complexity of the region of interest, it may 

15 be desirable to reduce to size of the dataset of the three dimensional volume prior to 
generating the skeleton. To this end, a multiresolution data reduction process, which is 
discussed in more detail below, can be used (step 341 5). 

The skeleton is a subset of the three dimensional volume. Preferably, 
the skeleton has the following attributes: 1) It preserves the homotopy of the tree; 2) it 

20 is 26-connected; 3) it is one voxel thick; 4) it approximates the central axes of the 
branches; and 5) it is relatively smooth. The degree of homotopy is somewhat 
application specific. For example, in generating a skeleton of the colon lumen, the 
skeleton will generally be a single path from end to end, despite the presence of 
numerous haustral folds which can be several voxels deep. However, in the 

25 cardiovascular system and pulmonary system, a small offshoot from the root which is 
several voxels deep can represent a legitimate branch in the system. 

Returning to Figure 34, in the volume of interest, a root voxel is 
identified in the volume (step 3420). In performing virtual endoscopy, this can be 
performed manually based on an understanding of the geometry of the structure being 

30 evaluated. 
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A distance map can then be generated to identify the branches in the 
tree and the distances between the endpoints of the branches and the root voxel (step 
3425). A presumption applied in this method is that there exists one unique endpoint 
on each branch which exhibits the longest distance to the root of the tree. Figure 35 is 
a schematic diagram illustrating a 3x3x3 cubic voxel arrangement which is referred to 
as a 26-connected voxel cubic distance plate. In the center of this arrangement is a 
seed voxel 3505, which is assigned a distance weight of zero. Around the seed voxel 
3505 are 26 connected neighbor voxels which are assigned distance weights based on 
the respective Euclidean distance between the respective neighbor voxel and the seed. 
In a cubic arrangement the Euclidian distance can assume a normalized value of 

'I -v/2, VS which is approximately equal to 1, 1.4 and 1.7. To simplify processing, 

the voxels can be assigned integer value weights of 10, 14, and 17 to approximate the 
relative Euclidian distances. 

Figure 36 is a pseudo-code representation of an algorithm for 
determining the distance map from a voxel in the volume to the root using the 
Euclidian weighted distances of the 26-connected cubic distance plate of Figure 35. 
From the generated distance map, branches are identified and the endpoints of the 
branches are determined (step 3430). 

Referring to Figure 36, the root of the volume is labeled with the integer 
value 0. A processing queue is then formed with the voxels in the volume. The 
voxels are then labeled in a first-in, first out manner by adding the Euclidian distances 
between the voxel at the top of the queue and the root voxel. This process is repeated 
until all of the voxels in the volume are assigned a value in the distance map. 

Because the labeling of voxels in the distance map will depend, in part, 
on the queuing order, the resulting distance map does not provide a unique solution. 
However, regardless of the queuing order, there is always at least one farthest point for 
each branch. In addition, for each voxel, other than the root voxel, there is always at 
least one 26-connected neighbor in the volume which has a shorter distance to the root. 
Thus, the endpoints of the branches are readily detectable by searching the distance 
m^ for local maximum distances (local maxima) (step 3430). The term local maxima 
is a relative term. In evaluating the volume for local maxima, the volume should be 
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partitioned into various subspaces which are appropriate to the object being evaluated. 
The expected feature size, branch length, branch diameter, etc, are generally considered 
in determining the subspace partitions. 

Once the endpoints of the branches are determined, the shortest path 
from the endpoint to the root voxel is determined (step 3435). The shortest paths from 
the endpoints to the root define the basic structure of the branches of the tree and 
approximate the centerline of the branches. This is referred to as the rough skeleton of 
the volume. The shortest paths are preferably generated from the branches at farthest 
end of the tree and begin from the end of those branches. From the most remote 
branch endpoint, the first voxel is selected and its 26-connected neighbors are analyzed 
to determine which voxel is in the minimal distance path from endpoint to root. This 
process is repeated until a selected voxel meets the root. This results in a one-voxel 
wide path from the farthest end to the root. Searching for the shortest path for other 
branches is sinoilar. However, for subsequent branches, the selection process can 
terminate when the current path reaches a previously assigned path (e.g., the path need 
not lead all the way to the root). The collection of all of the interconnected shortest 
paths is the resulting rough skeleton of the object. 

Depending on the application of the resulting rough skeleton, it may be 
desirable to refine the rough skeleton (step 3440). One step of refining the skeleton is 
to centralize the skeleton within the branches. Centralization preferably takes place 
branch by branch from the longest branch to the shortest. Starting with the longest 
branch, a uniform interval is selected, generally in the range of 4-8 voxels, along the 
branch. For each interval, the tangent direction of the voxel on the rough skeleton is 
calculated and a plane crossing the voxel perpendicular to the tangent direction is 
determined. A two dimensional area defined by the intersection of the plane and the 
volume is created and the center of this intersection can be computed using the known 
maximum disk technique. The centers of intersection can then be connected using a 
bi-cubic, B-spline interpolation or other curve fitting method. For the remaining 
branches, the endpoint which meets another branch or the root must first be adjusted to 
match the position of the previously centered skeleton branch. Then, centralization can 
proceed in the same manner as described for the longest branch. 
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Referring back to step 3415, when a large dataset is involved, it may be 
required, or at least desirable, to reduce the size of the dataset to speed up processing 
and reduce processing cost. Noting that the tree structure can be preserved within a 
range of scales, the large volume can be shrunk to a smaUer scale space for structure 
5 analysis. 

A shrinking method based on multiresolution analysis theory can be 
used. The input data is the stack of binary images of the same size which can be 
obtained from the segmentation results of the CT or MRI scan. The x-direction is 
taken along the slice image width, the y-direction is along the slice image height, and 

0 the z-direction is along the direction of slice by slice. The foreground voxels in the 
tree volume are set to value of 128 (maximum) and the background voxels are set to 
value 0 (minimum). A Daubechies' bi-orthogonal wavelet transform with all rational 
coefficients is employed. This one-dimensional (ID) discrete wavelet transformation 
(D WT) is first applied along to the x-direction row by row. From application of the 

5 • DWT only the lower frequency components are retained and packed. The computation 
is preferably implemented in floating points. Noting that the DWT is applied to the 
binary signal, there are two kinds of nonzero coefficients which result in the lower 
frequency component. The first is of value 128 and this kind of coefficients are located 
in the interior of the volume. The second is of a value not equal to 128 and these 

0 coefficients locate the boundary of the volume. 

The coefficients of the second kind are compared against a 
predetermined threshold value. If its absolute value is larger than a pre-set threshold 
Tl, the value of the coefficient is set to 128; otherwise, it is set to 0. This results in a 
stack of binary images with a row size of half of the original dataset. The same DWT 

5 is then applied to the resulting dataset along the y-direction column by column, where 
the similar thresholding is employed to the lower frequency components. The result is 
again a stack of binary images, but now with both half row and column size as 
compared to the original dataset. Finally, the DWT is applied to the last result along 
the z-direction and the lower frequency components are retained. This step completes 

0 the first level decomposition. 
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The resulting dataset of the first level decomposition is of half size in all 
three directions as compared to the original dataset. If the shrinking procedure stops at 
this level, the finial thresholding is applied. It revalues those coefiBcients of nonzero 
and non-128 value. If the absolute value of this kind of coefficient is larger than a pre- 

5 set threshold T2, it will be revalued as 128; otherwise, it is revalued as 0. If further 
shrinking is needed, the same thresholding algorithm is applied with the threshold Tl. 
Fxirther shrinking proceeds as previously described, but is applied to the dataset shrunk 
at the last previous level. The decomposition procedure can be recursively applied 
until the resulting volume meets the desired reduced data volume. In virtual 

10 endoscopy, the slice images are of 512X512 pixel size. The maximum decomposition 
level is usually three, resulting in a 64x64 reduced pixel size. 

The volume is isotropically shrank in all directions with the presented 
method. The two pre-set thresholds, Tl and T2, are used to control the degree of 
shrinking. If the volume is significantly over shrunk, connectivity may be lost in the 

15 reduced volume. If it is over shrunk to a leaser degree, two separate branches may 
merge into one branch in the reduced volume dataset. The larger the two threshold 
values, the thinner the reduced volume is. The range of those two thresholds is [0, r x 
128], where 0<r<l. Preferably, the range for virtual endoscopy is re (0.08, 0.28) for 
Tl and re (0.7, 0.98) for T2. The exact determination is dependant on the feature size 

20 of the particular application and is selected to achieve reduction while retaining the 
fidelity of the structure information in the shrunk volume. 

After shrinking the original volume, the tree branch searching procedure 
can be applied to the smaller volxune (steps 3420-3440). The resultant skeleton can be 
mapped back into the original scale space. When scaled to the original space, the 

25 image of the smaller scale skeleton no longer remain a connected path in the original 

* 

scale space. These voxels in the image act as control points for the fmal skeleton. The 
control points are centralized using the algorithm as described previously, and then, 
they are interpolated to form the final skeleton of the object. 

Computer Assisted Diagnosis 
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The virtual examination techniques described herein lend themselves to 
applications for the computer assisted diagnosis (CAD) of various conditions. For 
example, as described above, by examining the geometry of an organs tissue for local 
Gausian curvatures, regions with abnonnal geometry, such as polyps inside a colon 
S lumen, can automatically be identified. This technique can be generalized and used in 
conjunction with texture features to provide CAD functionality for a number of 
applications. 

For example, using the multi-scan imaging of the bladder described 
above, automated detection or tumors in the bladder wall can be performed. In this 

10 case, the degree of tumor invasion within the bladder wall is generally used to define 
the stage of bladder cancer. Using the multi-scan imaging and image segmentation 
techniques described above, the region of the bladder wall can be readily delineated. 
Regions of normal bladder tissue generally exhibit a substantially uniform texture 
feature. However, if a tumor is present in the region, the uniform texture feature will 

IS be interrupted. Thus, using texture analysis to evaluate the wall of the bladder, a 

region which may exhibit a tumor will present itself as a disturbance, or "noisy region" 
within the uniform texture. 

The texture of a region can be represented by a probabiUty distribution 
function (PDF) which characterizes the intensity correlation between voxels within a 

20 defmed range. A two-dimensional PDF can be used to represent a texture feature. 
Such a PDF characterizes the correlation between two closest voxels along all 
directions. To estimate the PDF, the intensities of any two closest neighbor voxels in a 
region of interest can be recorded as a sample vector for the region of interest (e.g., 
context). Using a number of such sample vectors, a cumulating distribution function 

25 (CDF) can be generated which estimates the PDF for that context. For each voxel, 

sample vectors within a range of its neighbor can also be used to generate a local CDF. 

A statistical test, such as a Kolmogorov-Smimov test, can be applied to 
the CDF to determine whether the CDF of the context and the local CDF are 
statistically equivalent, e.g., within a predefined confidence level. If so, the local 

30 texture feature around the current voxel is regarded as identical to the context. 
Otherwise, the current voxel exhibits a different texture feature from that of the 
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context and may be regarded as a potential abnormality, such as a tumor. The level of 
confidence used to determine whether a voxel is statistically equivalent to the context 
can be varied to increase or decrease the sensitivity of detection. 

In an alternative method of applying the PDF or CDF for texture 
anaiyis, each CDF or PDF can be regarded as a point in a functional linear space. The 
distance between two CDF's or PDF's in that space can be measured, such as in terms 
of the Skorohold metric. This distance provides«a measure of the degree of similarity 
of PDF's. For example, the distance between a local CDF and the context CDF can be 
calculated and the resulting distance can be compared to one or more distancing 
thresholds. If the distance is large, the local texture may be considered different from 
the context, which can indicate that such a voxel belongs to a region with a potential 
abnormality or tumor. Preferably, the distancing thresholds are determined based on 
evaluation of a statistically sufficient known data sets. 

The distance calculated above can be used v^th visualization techniques 
and volume rendering techniques, such as those described herein. For example, a 
feature volume dataset having a size comparable to the original dataset can be created. 
The intensity for each voxel in the new dataset can then be assigned based upon the 
distance between the local CDF and the CDF of the context. When this three 
dimensional volume dataset is viewed through volume rendering techniques, the 
regions which contain suspected tumors will exhibit a higher image intensity than the 
surrounding area. 

As was discussed above in connection with automatic detection of 
polyps, the surface of a lumen can be represented as a continuously second 
differentiable surface in three dimensional Euclidean space, such as by using a C-2 
smoothness surface model. In such a model, each voxel on the surface of the colon has 
an associated geometrical feature which has a Gauss curvature, referred to as Gauss 
curvature fields. For various organs, certain expected local features can be 
characterized by distinct curvature templates. For example, in the context of the colon, 
the expected local features include smooth curve surfaces, ring folds, convex hills from 
a smooth surface and plateaus from a smooth siurface. These last two local features 
may be indicative of a polyp or tumor. Accordingly, by searching the Gauss curvature 
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fields for specific predetermined local feature templates, polyps, tumors and other 
abnormalities of interest can be automatically detected. This use of surface geometry 
to perform computer assisted diagnosis can be used alone, or in conjunction with the 
texture-based CAD techniques described above. 
5 As an alternative to rendering and viewing the entire organ or region of 

* 

interest, the surface under observation can be partitioned into small areas, or patches, 
which are defined by the local curvature templates. Each patch should contain voxels 
which have a common geometry feature, or curvature template. A single viewing point 
is then determined for the patch which allows all voxels of the patch to be observed. 

0 The patches are then assigned a priority score indicating the probability that the patch 
represents a polyp or other abnormality. The patches can then be observed 
individually, in priority order, rather than requiring the operator to navigate the entire 
organ volume to search out suspect areas. Of course, a preferred diagnostic system 
includes the ability to toggle between views such that an operator can readily change 

5 firom viewing a patch to viewing the patch in the context of the organ. Alternatively, 
these two views can be presented simultaneously. Again, the texture based approaches 
can be used to supplement this approach. By mapping the results of texture analysis 
onto the patches being observed, the texture information can also be observed and used 
in diagnoses. 

0 The foregoing merely illustrates the principles of the present imaging 

and examination systems and methods. It will thus be appreciated that those skilled in 
the art will be able to devise numerous systems, apparatus and methods which, 
although not explicitly shown or described herein, embody the principles of the 
invention and are thus within the spirit and scope of the invention as defined by its 

5 claims. 

For example, the methods and systems described herein could be 
applied to virtually examine an animal, fish or inanimate object. Besides the stated 
uses in the medical field, applications of the technique could be used to detect the 
contents of sealed objects which cannot be opened. The techniques can also be used 
0 inside an architectural structure such as a building or cavern and enable the operator to 
navigate through the structure. 
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CLAIMS 

1 . A method for performing virtual examination of an object comprising: 

performing at least one scan of an object with the object distended by 
the presence of a contrast agent; 
5 performing at least one scan of the object with the object relieved of the 

contrast agent; ^ 

converting the scans to* corresponding volume datasets comprising a 

plurality of voxels; 

performing image segmentation to classify the voxels of each scan uito 
10 a plurality of categories; 

registering the volume datasets of each scan to a common coordinate 

system; 

displaying at least two of the volume datasets in a substantially 
simultaneous manner; and 
1 5 performing virtual navigation operations in one of the volume datasets 

and having the corresponding navigation operations take place in at least one other 
volume dataset. 

2. The method for performing virtual examination according to claim 1 , wherein 
the at least one scan of the distended object includes a transverse scan and a coronal 

20 scan of the object 

3. The method for performing virtual examination according to claim 2, wherein 
the at least one scan of the relieved object includes a transverse scan and a coronal scan 
of the object. 

4. The method for performing virtual examination according to claim 3, wherein 
25 the object is a bladder. 
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5. The method of perfonning virtual examination according to claim 4, wherein 
the scans are computed tomography scans. 

6. The method of performing virtual examination according to claim 4, wherein 
the scans are ultrasound imaging scans. 

5 7. The method of perfonning virtual examination according to claim 4, wherein 
the scans are magnetic resonance imaging scans. 

* 

8 The method of performing virtual exammation according to claim 7, wherein 
the contrast agent is urine. 

9. The method for performing virtual examination according to claim 1, wherein 
10 the at least one scan of the relieved object includes a transverse scan and a coronal scan 

of the object. 

10. The method for performing virtual examination according to claim 1, wherein 
the object is a bladder. 

1 1 . The method of performing virtual examination according to claim 1 0, wherein 
15 the scans are computed tomography scans. 

12. The method of performing virtual examination according to claim 10, wherein 
the scans are ultrasound imaging scans. 

13. The method of perfonning virtual examination according to claim 10, wherein 

* • 

the scans are magnetic resonance imaging scans 

20 14. The method of performing virtual examination according to claim 1 3, wherein 
the contrast agent is urine. 
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15. The method of performing virtual examination according to claim 1 , further 
comprising evaluating the at least one scan with the object distended and the at least 
one scan with the object relieved to identify regions where contrast is more visible in 

■ 

one of said scans and evaluating the scan with more contrast in a region of interest to 
5 determine physiological characteristics of the object. 

16. The method of performing virtual examination according to claim 15, wherein 
said step of image segmentation includes classifying voxels based on local intensity 
vectors of the voxels. 

17. The method of performing virtual examination according to claim 16, wherein 
10 the step of image segmentation further includes using a region growing algorithm to 

identify regions of the object based on the classified voxels. 

1 8. The method of performing virtual examination according to claim 1 , further 
comprising partitioning the volume image datasets into a plurality of regions related to 
the coordinate system. 

15 19. The method of performing virtual examination according to claim 1 8, wherein 
the plurality of regions include eight regions defined in a three dimensional coordinate 
system. 

20. - A method for performing virtual examination of an object comprising: 
20 performing an imaging scan of the object to acquire image scan data; 

converting the acquired image scan data to a plurality of voxels; 
interpolating between the voxels to generate an expanded dataset; 
performing image segmentation to classify the voxels into a plurality of 

categories; 

25 extracting a volume of the object interior from the expanded dataset; 

generating a reduced resolution dataset from the expanded dataset; 
storing the expanded dataset in a tree data structure; 
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rendering images for the expanded dataset and reduced resolution 

dataset; and 

selecting at least one of the reduced resolution dataset or expanded 
dataset renderings for display. 



5 21. The method for performing virtual examination of an object of claim 20, 
wherein the selecting step comprises: 

selecting the reduced resolution dataset during image interaction; and 
selecting the expanded dataset rendering if no image interaction has 
occurred in a predetermined time period. 



1 0 22. The method for performing virtual examination of an object of claim 20, 
wherein the imaging scan is a computed tomography scan. 



23. The method for performing virtual examination of an object of claim 20, 
wherein the imaging scan is a magnetic resonance imaging scan. 

24. The method for performing virtual examination of an object of claim 20, 
1 5 wherein the imaging scan is an ultrasound imaging scan. 

25. The method for performing virtual examination of an object of claim 20, 
wherein the object is the larynx, 

26. The method for performing virtual examination of an object of claim 20, 
wherein the tree structure is a binary space partition tree structure. 

20 27. A method of performing virtual angiography comprising: 

acquiring imaging scan data including at least a portion of the aorta; 
converting the imaging scan data to a volume representation including 
plurality of voxels; 
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segmenting the volume representation to classify the voxels into one of 
a plurality of categories; 

analyzing the segmented volume representation to identify voxels 
indicative of at least a portion of an aneurysm in the aortic wall; and 
5 generating at least one closing surface around the voxels indicative of at 

least a portion of an aneurysm to estimate the contour of the aneurysm. 

28. The method of performing virtual angiography of claim 27, wherein the 
imaging scan is a computed tomography scan. 

29. The method of performing virtual angiography of claim 27, wherein the 
1 0 imaging scan is a magnetic resonance imaging scan. 

30. The method of performing virtual angiography of claim 27, wherein the 
segmenting operation classifies voxels in at least the categories of blood, tissue, and 
calcium deposits. 

31. The method of performing virtual angiography of claim 27, further comprising 
1 5 estimating the volume of the aneurysm using the generated closing surfaces. 

» • 

32. The method of performing virtual angiography of claim 27, further comprising 
generating a navigation path through the aortic lumen. 

33. The method of performing virtual angiography of claim 32, further comprising 
20 estimating the length of the aneurysm based on the navigation path. 

34. A method of performing virtual endoscopy of a blood vessel comprising: 

acquiring imaging scan data including at least a portion of the vessel; 
converting the imaging scan data to a volume representation including a 
plurality of voxels; 
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segmenting the volume representation to classify the voxels into one of 
a plurality of categories including the categories of blood, tissue, and calcium deposits; 
and 

generating a navigation path through the vessel. 

5 35. The method of performing virtual endoscopy of claim 34, wherein the vessel is 
a carotid artery. 

36. The method of performing virtual endoscopy of claim 34, further comprising 
the step of determining the diameter of the carotid artery along the navigation path to 
identify regions of narrowing. 

10 37. The method of performing virtual angiography of claim 34, wherein the 
imaging scan is a computed tomography scan. 

38. The method of performing virtual angiography of claim 34, wherein the 
imaging scan is a magnetic resonance imaging scan. 

39. A method of determining the characteristics of a stent graft using virtual 

1 5 angioscopy, comprising: 

acquiring imaging scan data including at least a portion of the aorta; 
converting the imaging scan data to a volume representation including a 

plurality of voxels; 

segmenting the volume representation to classify the voxels into one of 

20 a plurality of categories; 

analyzing the segmented volume representation to identify voxels 
indicative of at least a portion of an aneurysm in the aortic wall; 

generating at least one closing surface around the voxels indicative of at 
least a portion of an aneurysm to estimate the contour of the aneurysm; 
25 identifying the location of the endpoints of the aneurysm contour; 
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calculating the length between the endpoints of the aneurysm contour to 
determine the length of the stent graft; and 

calculating the diameter of the aortic lumen at the endpoints of the 
aneurysm contour to determine the required outside diameters of the stent graft. 

5 40. The method of determining the characteristics of a stent graft of claim 39, 
fijrther comprising determining the angle of interface of the aneurysm and normal 
aortic lumen to determine an angular direction of a corresponding end of the stent 
graft. 

41 . The method of determining the characteristics of a stent graft of claim 39, 
10 ftirther comprising locating arterial branches proximate the aneurysm to determine a 

maximum length of the stent graft. 

42. The method of determining the characteristics of a stent graft of claim 4 1 , 
wherein the arterial branches proximate the aneurysm include at least one of the renal 
and femoral arterial branches. 

1 5 43. The method of determining the characteristics of a stent graft of claim 39, 

ftuther comprising conducting a virtual biopsy of the aortic region proximate the ends 
of the aneurysm to determine the nature of the tissue at the anticipated graft interface 
locations. 

44. A method of defining a skeleton for a three dimensional image representation 
20 of a hollow object formed with a plurality of voxels comprising: 

identifying a root voxel within the hollow object; 

generating a distance map for all voxels within the hollow object, the distance 
map being formed using a 26-coimected cubic plate of neighboring voxels having 
Euclidian weighted distances; 
25 identifying voxels having a local maxima in the distance map as endpoints of 

branches in the hollow object; and 
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for each local maxima voxel, determining a shortest connected path to one of 
the root voxel or a previously defined shortest path. 

45. The method of defining a skeleton for a three dimensional image representation 
of claim 44 further comprising performing multi-resolution data reduction to the three 

5 dimensional image representation to generate a reduced data set for the generating and 
identifying operations. 

46. The method of defining a skeleton for a three dimensional image 
representation of claim 44 further comprising centralizing the shortest paths within the 
respective branches of the object 

1 0 47. The method of defining a skeleton for a three dimensional image 

representation of claim 44, wherein the object mcludes at least one blood vessel. 

48. The method of defining a skeleton for a three dimensional image representation 
of claim 44, wherein the object includes the airways of a Ixmg. 

49. The method of defming a skeleton for a three dimensional image representation 
1 5 of claim 44, wherein the object includes the bladder. 

50. The method of defining a skeleton for a three dimensional image representation 
of claim 44, wherein the object includes the spinal cord of a vertebrate animal. 

51. A method of performing computed assisted diagnosis of a region of interest, 
comprising: 

20 acquiring imaging scan data including at least a portion of the region of 

interest; 
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converting the imaging scan data to a volume representation including a 
plurality of voxels, at least a portion of the voxels representing a surface of the region 
of interest; and 

analyzing said portion of voxels representing a surface for at least one 
5 of a geometric feature and a textural feature indicative of an abnormality. 



52. The method of performing computed assisted diagnosis according to claim S 1 , 
wherein the textural feature is included in a probability density function characterizing 
a correlation between two voxels of the portion of voxels. 



53 The method of performing computed assisted diagnosis according to claim 52, 
1 0 wherein the two voxels are adj acent voxels. 



54. The method of performing computer assisted diagnosis according to claim 52, 
wherein intensities of said portion of voxels are used to generate an estimate of the 
probability density function. 



55. The method of performing computer assisted diagnosis according to claim 54, 
1 5 wherein a plurality of voxel intensities are used to generate a cumulating distribution 
function of the region of interest and a local cumulating distribution function, and 
wherein the local cimiulating distribution function is compared against the context 
cumulating distribution function to identify regions of abnormality. 



56. The method of performing computer assisted diagnosis according to claim 55, 
20 wherein a distance is determined between said local cumulating distribution function 
and said context cumulating distribution function, the distance providing a measure of 
abnormality. 



57. The method of performing computer assisted diagnosis according to claim 56, , 
wherein the distance is used to assign intensity values to the voxels representing a 
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surface of the region of interest and wherein said method further comprises displaying 
said voxels such that variations in intensity represent regions of abnormality. 



58. The method of performing computer assisted diagnosis according to claim 57, 
wherein the region of interest includes the colon and. wherein the abnormality includes 

5 polyps. 

59. The method of perfonniiig computer assisted diagnosis according to claim 51, 
wherein the region of interest includes the aorta and wherein the abnormality includes 
abdominal aortic aneurysms. 

m 

60. The method of performing computer assisted diagnosis according to claim 5 1 
10 wherein the surface is represented as a second diflferentiable surface where each surface 

volume unit has an associate Gauss curvature and wherein said Gauss curvatures 
combine to form said geometric features. 

61. The method of performing computer assisted diagnosis according to claim 59 
wherein a plurality of predetermined geometrical feature templates are defined and 

1 5 wherein the geometric features of said surface are compared to said templates to 
determine a geometric feature classification. 



wo 02/29764 



1126 



PCT/USOl/30704 



101 



r 



1 



103 



104 



105 



107 



( START ) 



PREPARE ORGAN IF 
NECESSARY 



CONVERT SCAN DATA 
INTO VOLUME ELEMENTS 



I 



DEFINE PORTION OF 
ORGAN TO EXAMINE 



I 




PERFORM AUTOMATIC OR 
GUIDED NAVIGATION OF ORGAN 



109 



I 



DISPLAY ORGAN 
i 

( END J 

FIG. 1 




wo 02/29764 



2/26 



PCT/USOl/30704 




# 



wo 02^764 PCT/USOl/30704 

3/26 



401 




FIG. 6 FIG. 7 




FIG. 8 




wo 02/29764 



PCTAJSOl/30704 



4/26 



901 



903 



905 



911 



913 



915 



( START ) 



IDENTIFY CELL 
CONTAINING CAMERA 



1 



BUILD STAB TREE FOR 
POTENTIALLY VISIBLE CELLS 



I 



STORE INTERSECTION OF 
ADJOINING CELLS AT EDGES 
OF STAB TREE 




909 



COLLAPSE 
TWO NODES 



INIALIZE Z-BUFFER 



I 



TRAVERSE NODES IN Z-BUFFER 



I 



build image of 
visible cells 

(end ) 



FIG. 9 




wo 02/29764 



PCTAJSOl/30704 



5/26 



cells 



advancing cross -section 



portals 




finish point 



start point 



center-line 



FIG. 10 



1103 




FIG. 11(a) 



wo 02/29764 



6/26 



PCT/USOl/30704 



1109 




INTERSECT 
{A/B,B/C} 



1112 



A/B 




INTERSECT 
{A/B, B/D} 



INTERSECT 
{A/B, B/C, C/B} 




•1114 




1109 



A/B 




INTERSECT 
{A/B, B/D} 



1114 




FIG. 11(b) 



FIG. 11(c) 



1261 
1257 



1203 



1253 



1255 



1201 1267 



1269 




Km 6 



1251 1269 1203 



1265 



1263 



FIG. 12(a) 



wo 02/29764 



PCTAJSOl/30704 



7/26 



1211 




I/F 



1213 




INTERSECT, 
{I/F,F/By 



INTERSECT 
{l/F, F/E} 



1215 





1219 



INTERSECT 
{l/F, F/B, B/A} 



1217 




FIG. 12(b) 



INTERSECT; 
{I/F, F/B}^ 

1215 





INTERSECT 
{I/F, F/E} 



1219 



1217 




I/F 



1215 




INTERSECT 
{I/F. F/B} 



INTERSEa 
{I/F, F/E} 



1219 



INTERSECT 
(l/F, F/B, B/a} 



1217 




RENDERED NODES {l} 
SKIPPED N0DE{} 

FIG. 12(c) 



INTERSECT 
{I/F, F/E} 



1219 



INTERSECT 
{I/F, F/B, B/A> 

RENDERED NODES: {I, F} 
SKIPPED NODES: { } 

FIG. 12(d) 




RENDERED NODES: {l,F} 
SKIPPED NODES:{A,B} 

FIG. 12(e) 



'7 



1301 



1309 
1307 / 





wo 02/29764 



8/26 



1405< 




1401 



1402 



1407 



X 



1406 



/ 



CPU 



-1409 



MEMORY 



•1411 



•1413 



1427 



■1417 



1416< 



1415 



CPU 



-1423 



MEMORY 



■1425 



PCT/USOl/30704 



1414 




■1421 



1431 



(1111 



1 



9 



c 




1419 



r 



Fin 14. 



wo 02/29764 



9/26 




START 




REGIONAL 
PREPARATION FOR 
IMAGING 



PCT/USOl/30704 



■1510 



ig. 15 



ACQUIRE IMAGE 
DATA 



I 



1520 



CONVERT IMAGE 
DATA TO VOXELS 



^1530 



GROUP VOXELS BY 
INTENSITY 
THRESHOLD AS 
VOXEL CLUSTERS 



IDENTIFY TARGET 
CLUSTER FOR 
REGION OF 
INTEREST 



PERFORM FEATURE 
VECTOR ANALYSIS 
ON TARGETED 
CLUSTER<S) 



CLASSIFY VOXELS 
WITHIN CLUSTER 



PERFORM HIGH 
LEVEL FEATURE 
EXTRACTION 



1540 



1550 



1560 



1570 



_^580 



,. t 

END 



wo 02/29764 



10/26 



PCT/USOl/30704 




The fixed voxe 



The second nearest slice 
The fixed slice 




wo 02/29764 



11/26 



PCT/USOl/30704 




f,6r. ISC 



wo 02/29764 



12/26 



PCT/USOl/30704 



Fig. 20 



Start 





• 


TEXTURE 
SEGMENTATION 






TEXTURE 
ANALYSIS 







TEXTURE 
MODELING 






TEXTURE 
MATCHING 







TEXTURE 
SYNTHESIS 



2010 



'2020 



■2030 



2040 



2050 



END 



wo 02/29764 



FIG. 21 



13/26 



CAST RAY FROM 
THE VIEWPOINT 
THROUGH IMAGE 
PIXELS 



PCT/USOl/30704 



SET FIRST 
SAMPLING POINT AS 

THE CURRENT 
IMAGE PIXEL ALONG 
THE RAY 





CHECK DISTANCE 
FROM CURREhfT 
SAMPLING POINT TO 
NEAREST COLON 
WALL 


4 



















JUMP TO NEW 
POINT ALONG RAY 
WITH DISANCE (D) 



Yes- 



DISTANCE (D)> 
SAMPLING INT. (1)7. 




PERFORM REGULAR 
SAMPLING AT THIS 
POINT 



r 



GOTO NEW 




SAMPLING POINT 




ALONG RAY AT 




DISTANCE i 





wo 02/29764 



14/26 



PCT/USOl/30704 



FIG. 22 



START 



VOLUME SHRINKING 
FROM COLON WALL 



GENERATE MIN. 

DISTANCE PATH 
BETWEEN 
ENDPOINTS 



EXTRACT CONTROL 
POINTS 



CENTER CONTROL 
POINTS WITHIN 
LUMEN 



INTERPOLATE 

POINTS TO 
GENERATE FLY- 
PATH THROUGH 
LUMEN 



STOP 



wo 02/29764 



15/26 



PCT/USOl/30704 



START 



FIG. 23 



REPRESENT 
VOLUME AS A 
STACK OF BINARY 
IIMAGES 



1- 



V 



APPLY DISCRETE 

WAVELET 
TRANSFORMATION 



"Vi 



APPLY THRESHOLD 
TO SUB-DATASETS 





wo 02/29764 



PCT/USOl/30704 



16/26 




wo 02/29764 



PCT/USOl/30704 



17/26 



START 



FIG. 26 



YES 




SEGMENT COLON 
LUMEN . 



2510 



SELECT POINT 
WITHIN EACH 
. SEGMENT 



2520 




CENTER SELECTED 

POINTS WITH 
RESPECT TO COLON 
WALL 



2530 



SEQUENTIALLY 
CONNECT 
CENTERED POINTS 



2540 



-I ^ 2550 

END 



wo 02/29764 



PCT/USOl/30704 



18/26 



IMAGE DATA FROM 
1AAA6ING SCANNER 




SCANNER 
INTERFACE 
BOARD 



/-A 
W 






VOLUME 


VOLUME 


RENDERING 


RENDERING 


MEMORY 


BOARD 





M 
W 



GRAPHICS 
BOARD 



W 



MAiN MEMORY 



I/O DEVICES 



wo 02/29764 



19/26 



PCT/USOl/30704 



CONVERT IMAGING 
SCANNER DATA TO 
VOLUME ELEMENT 
REPRESENTATION 
(MAIN MEMORY) 



^1/ 




FOR EACH SLICE, 
PERFORM VOLUME 
RENDERING AND 
TEXTURE SYNTHESIS 
(VOLUME RENDERING 
MEMORY) 



STORE RENDERED 
SLICES IN SEQUENTIAL 
BUFFER (BACK TO 
FRONT) 
(MAIN MEMORY) 



DISPLAY CONTENTS OF 
BUFFER 



wo 02/29764 



20/26 



PCT/USOl/30704 



2805 




2810 



2815- 



2820 



2825 



2830 



Pre-lmage 
SCAN Protocol 



r 





Scan 1 




Bladder full, 




transverse 




scan 



Scan 2 
Bladder full, 
Coronal scan 



Scans 
Bladder 
Relieved, 
transverse 
scan 




I 



Scan 4 
Bladder 
Relieved, 
Coronal scan 



IMAGE 
SEGMENTATION 



2835 



REGISTER 
COORDINATE 
FRAMES OF 

OCTANTS 



2840 








PARTITION 




DATA SETS 




(OCTANTS) 


2845 








VIEW MULTI- 




SCAN IMAGE 




wo 02/29764 



PCTAJSOl/30704 



21/26 






wo 02/29764 



22/26 



PCT/USOl/30704 



SCAN REGION OF 
LARYNX 



I 



3105 



F-1 pi 



IMAGE 
SEGMENTATION 

(GENERATE 
EXPANDED DATA 
SET AND CLASSIFY 
VOXELS) 



3110 




1 



EXTRACT VOLUME 
OF LARYNX 
VOLUME 



r 



3115 



GENERATE 
REDUCED 
RESOLUTION 
DATASET 



3120 



I 



EXTRACT 
POLYGONS 
REPRESENTING 

SURFACES/ 
VOLUMES FROM 
BOTH DATASETS 



T 



3125 



ORGANIZE 
ENLARGED 
DATASET IN A BSP- 



3130 



DISPLAY LOW- 
RESOLUTION OATSET 



4 — Yes 



15 




3135 



DISPLAY HIGH- 
RESOLUTION 
DATASET 



3140 



wo 02/29764 



23/26 



PCT/USOl/30704 



'G. 32 




3205 



3210 




IDENTIFY' 

PORTIONS OF 
ANEURISM 



GENERATE CLOSING 
SURFACE 



3215 



3220 



3225 



ESTIMATE VOLUME 
OF ANEURISM 



wo 02/29764 



24/26 



PCTAISOl/30704 



/V^ 33A 




wo 02/29764 



25/26 



PCT/USOl/30704 



Acquire image Scan 



3405 



Generate 3-D Volume 
of Object 



3410 



- 



Identify Root 



3415 



Multiresolution Data 
Reduction 



3420 



X 



Generate Distance 
Map to the Root 



3425 



Determine Ends of 
Branches 



3430 



Search for Shortest 
Paths 



3435 



3440 



Refine Rough Skeleton 



wo 02/29764 



26/26 



PCT/USOl/30704 





1) Label root voxel with integer 0; 

2) Construct a queue and line up the root in the queue; 

3) If( There is at least one voxel in the queue ) 

Seive the voxel x on the top of the queue: 
For( each ofx's26-connected neighbor voxel y) { 
If( y in the volume and has not been labeled yet ) { 
Line up the y in the queue; 

/* label the voxel y */ 
Set dist = 999999 ; 

For( each of y's 26<onnected neighbor voxel z ) { 

If( z in the volume and has been already labeled with an integer of ) ( 

where cHy,z)is 10, 14 or 17 if the Euclidean disunce between y and z is 
1.^2^ or V3 . respectively; 

) 

1t(.dist>d^) { 

label y with integer dist ; 
dist = ; 

) 

) 

) 

) 

X leaves the queue; 

} 

Else { 

end of calculating the distance m^. 

) 
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SYSTEM AND METHOD FOR PERFORMING A 
THREE-DIMENSIONAL VIRTUAL EXAMINATION OF OBJECTS, SUCH 

AS INTERNAL ORGANS 

SPECraCATION 

TECHNICAL FIELD 
The present invention relates to a system and method for performing a 
volume based three-dimensional virtual examination, and more particularly relates to 
a system which offers enhanced visualization and navigation properties. 

BACKGROUND OF THE INVENTION 
Colon cancer continues to be a major cause of death throughout the 
world. Early detection of cancerous growths, which hi the human colon initially 
manifest themselves as polyps, can greatly improve a patient's chance of recovery. 
Presently, there are two conventional ways of detecting polyps or other masses in the 
colon of a patient. The first method is a colonoscopy procedure, which uses a flexible 
fiber-optic tube called a colonoscope to visually examine the colon by way of physical 
rectal entry with the scope. The doctor can manipulate the tube to search for any 
abnormal growths in the colon. The colonoscopy, although reliable, is both relatively 
costly in money and time, and is an invasive, uncomfortable painful procedure for the 
patient. 

The second detection technique is the use of a barium enema and two- 
dimensional X-ray imaging of the colon. The barium enema is used to coat the colon 
with barium, and a two-dimensional X-ray image is taken to capture an inciage of the 
colon. However, barium enemas may not always provide a view of the entire colon, 
require extensive pretreatment and patient manipulation, is often operator-dependent 
when performing the operation, exposes the patient to excessive radiation and can be 
less sensitive than a colonoscopy. Due to deficiencies in the conventional practices 
described above, a more reliable, less intrusive and less expensive way to examine the 



wo 02/029764 



PCT/USOl/30704 



2 

colon for polyps is desirable. A method to examine other hmnan organs^ such as the 
lungs, for masses in a reliable, cost effective way and with less patient discomfort is 
also desirable. 

Another leading cause of cancer deaths in the United States is bladder 

5 cancer. In 1995, there were 50,000 new cases of bladder cancer reported and 1 1,000 
deaths were reported as a result of this disease. The most common test for bladder 
cancer is the use of a urine "dipstick" or conventional urinalysis. However, such tests 
are generally only effective at detecting bladder cancer in its later developed stages 
and does not provide any information regarding the size or location of a cancerous 

0 growth. Cystoscopy, the main method of investigating bladder abnoimalities at 
present, provides accurate results and can provide information regarding the relative 
size and location of any abnormalities. However, cystoscopy is an invasive procedure 
which offers a physician a limited field of view and lacks an objective indication of 
size. In addition, cystoscopy is contra-indicated for those patients who have severe 

5 urethral strictures or active vesical bleeding. Thus, it is desirable to develop 

alternative procedures for screening patients for bladder cancer, especially at early 
stages of cancer development. 

Two-dimensional ("2D") visualization of human organs employing 
currently available medical inaiaging devices, such as computed tomography and MRI 

0 (magnetic resonance imaging), has been widely used for patient diagnosis. Three- 
dimensional images can be formed by stackiag and interpolating between two- 
dimensional pictures produced firom the scanning machines. Imaging an organ and 
visualizing its voliune in three-dimensional space would be beneficial due to its lack 
of physical intrusion and the ease of data manipulation. However, the exploration of 

5 the three-dimensional volume image must be properly performed in order to fully 
exploit the advantages of virtually viewing an organ jfrom the inside. 

When viewing the three dimensional ("3D") volume virtual image of 
an environment, a functional model must be used to explore the virtual space. One 
possible model is a virtual camera which can be used as a point of reference for the 

0 viewer to explore the virtual space. Camera control in the context of navigation 
within a general 3D virtual enviroEment has been previously studied. There are two 
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conventional types of camera control offered for navigation of virtual space. The first 
gives the operator full control of the camera which allows the operator to manipulate 
the camera in different positions and orientations to achieve the view desired. The 
operator will ia effect pilot the camera. This allows the operator to explore a 
5 particular section of interest while ignoring other sections. However, complete 
control of a camera in a large domain would be tedious and tiring, and an operator 
might not view all the important features between the start and finishing point of the 
exploration. 

The second technique of camera control is a planned navigation 
1 0 method, which assigns the camera a predetermined path to take and which caimot be 
changed by the operator. This is akin to having an engaged "autopilot". This allows 
the operator to concentrate on the virtual space being viewed, and not have to worry 
about steering into walls of the environment being examined. However, this second 
^ technique does not give the viewer the flexibility to alter the course or investigate an 
1 S interesting area viewed along the flight path. 

0 

It would be desirable to use a combination of the two navigation 
techniques described above to realize the advantages of both techniques while 
minimizing their respective drawbacks. It would be desirable to apply a flexible 
navigation technique to the examination of human or animal organs which are 

20 represented in virtual 3D space in order to perform a non-intrusive painless thorough 
examination. The desired navigation technique would fiirther allow for a complete 
exairdnation of a virtual organ in 3D space by an operator allowing flexibility while 
ensuring a smooth path and complete examination through and around the organ. It 
would be additionally desirable to be able to display the exploration of the organ in a 

25 real time settuig by using a technique which minimizes the computations necessary 
for viewing the organ. The desired technique should also be equally applicable to 
exploring any virtual object. 

It is another object of the invention to assign opacity coefiBcients to 
each volume element in the representation in order to make particular volume 

30 elements transparent or translucent to varying degrees in order to customize the 
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visualization of the portion of the object being viewed. A section of the object can 
also be composited using the opacity coefiScients. 

SUMMARY OF THE INVEimO 
The invention generates a three-dimensional visualization image of an 

5 object such as a human organ using volume visualization techniques and explores the 
virtual image using a guided navigation system which allows the operator to travel 
along a predefined flight path and to adjust both the position and viewing angle to a 
particular portion of interest in the image away from the predefined path in order to 
identify polyps, cysts or other abnormal features in the organ. 

10 A method for performing virtual examination of an obj ect includes 

performing at least one imaging scan of an object with the object distended by the 
presence of a contrast agent In addition, at least one imaging scan of the object is 
acquired with the object relieved of the contrast agent The scans are converted to . 
corresponding volume datasets formed with a plurality of voxels. Image segmentation 

15 is then performed to classify the voxels of each scan into a plurality of categories. 
The volume datasets of each scan are registered to a common coordinate system. A 
displaying operation can then be performed where corresponding images at least two 
of the volume datasets are substantially simultaneously displayed. Virtual navigation 
operations performed m one of the volume datasets results in having the 

20 corresponding navigation operations take place in at least one other volume dataset. 

Preferably, the at least one scan of the distended object includes a 
transverse scan and a coronal scan of the object Similarly, it is preferable that the at 
least one scan of the relieved object includes a transverse scan and a coronal scan of 
the object. This procedure is particularly well suited for performing virtual 

25 cystoscopy, where the object is the bladder. In this case, the scan generally takes the 
form of a magnetic resonance imaging scan and the contrast agent can be urine. 

Another method in accordance with the present invention is for 
perforniing virtual exanunation of an object. In this method, an imaging scan of the 
object is performed to acquire image scan data. The acquired image scan data is 

30 converted to a plurality of volume units, or voxels. By interpolating between the 
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voxels^ an expanded dataset is generated. Image segmentation can then be performed 
to classify the voxels into a plurality of categories. A volume of the .object interior is 
extracted from the expanded dataset, such as by usrag a region growing algorithm 
from a seed voxel within the object lumen. A reduced resolution dataset is then 
5 generated from the expanded dataset. To efficiently store and recall the data from the 
expanded data set^ this data is stored in a tree data structure. Images can then be 
rendered for both the expanded dataset and reduced resolution dataset. One of these 
images is then selected for viewing. Generally, the reduced resolution dataset is 
selected during navigation or image interaction whereas the expanded dataset is 

a 
I 

1 0 selected for high resolution, static display. 

A method of performing virtual angiography is also provided. In this 
method, imaging scan data is acquired of at least a portion of the aorta. The imaging 
scan data is converted to a volume representation including a plurality of voxels. The 
volume representation is segmented to classify the voxels into one of a plurality of 

1 5 categories. The segmented volume representation is then analyzed to identify voxels 
indicative of at least a portion of an aneiuysm in the aortic wall. From the portions of 
the aneurysm which are identified, at least one closing surface is generated around the 
voxels indicative of at least a portion of an aneurysm. The closmg surface provides an 
estimate of the contour of the aneurysm. A navigation path can be established 

20 through the aortic lumen and characteristics of the aneurysm, such as length, diameter, 
volume and composition, can be determined. 

The method of performing virtual angiography can be used to detect 
and monitor the progression of aneurysms and can also be used in determining the 
characteristics needed to place a stent graft. 

25 Also provided is a method of defining a skeleton for a three 

dimensional image representation of a hollow object formed with a plurality of 
voxels. A root voxel is first identified within the hollow object. A distance map is 
then generated for all voxels within the hollow object The distance map is formed 
using a 26-connected cubic plate having Euclidian wei^ted distances. Those voxels 

30 haviag a local maxima in the distance map are identified as endpoints of branches in 
the hollow object. For each local maxima voxel, a shortest connected path to either . 
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the root voxel or a previously defined shortest path, is determined. The collection of 
shortest paths is the rough skeleton of the object. This technique is particularly well 
suited for multibxanch structures such as the respiratory system and cardio vascular 
system. 

5 

BRIEF DESCRIPTIQN OF THE DRAWINGS 
Further objects, features and advantages of the invention will become 
apparent firom the following detailed description taken in conjunction with the 
accompanying figures showing a preferred embodiment of the invention, on which: 

Figure 1 is a flow chart of the steps for performing a virtual 
examination of an object, specifically a colon, in accordance with the invention; 

Figure 2 is an illustration of a "submarine" camera model which 
performs guided navigation in the virtual organ; 

* 

Figure 3 is an illustration of a pendulum used to model pitch and roll of 
the "submarine" camem; 

Figure 4 is a diagram illustrating a two dimensional cross-section of a 
volumetric colon which identifies two blocking walls; 

Figure 5 is a diagram illustrating a two dimensional cross-section of a 
volumetric colon upon which start and finish volume elements are selected; 

Figure 6 is a diagram illustrating a two dimensional cross-section of a 
volumetric colon which shows a discrete sub-volume enclosed by the blocking walls 
and the colon surface; 

Figure 7 is a diagram illustrating a two dimensional cross-section of a 
volumetric colon which has multiple layers peeled away; 

Figure 8 is a diagram illustrating a two dimensional cross-section of a 
volumetric colon which contains the remaining flight path; 

Figure 9 is a flow chart of the steps of generating a volume 
visualization of the scanned organ; 

Figure 10 is an illustration of a virtual colon which has been sub- 
divided into cells; 
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Figure 1 1 A is a graphical depiction of an organ which is being virtually 

examined; 

Figure 11 B is a graphical depiction of a stab tree generated v^en 
depicting the organ in Fig. 1 1 A; 
5 Figure 1 IC is a further graphical depiction of a stab tree generated 

while depicting the organ isx Fig. 1 1 A; 

Figure 12A is a graphical depiction of a scene to be rendered with 
objects within certain cells of the scene; 

Figure 12B is a graphical depiction of a stab tree generated while 
0 . depicting the scene in Fig. 12 A; 

Figures 12C-12E are further graphical depictions of stab trees 
generated while depicting the image in Fig. 12A; 

Figure 13 is a two dimensional representation of a virtual colon 
containing a polyp whose layers can be removed; 
S Figure 14 is a diagram of a system used to perform a vutual 

examination of a human organ in accordance with the invention; 

Figure 15 is a flow chart depicting an improved image segmentation 

method; 

Figure 16 is a graph of voxel intensity versus frequency of a typical 
:0 abdommal CT data set; 

Figure 17 is a perspective view diagram of an intensity vector structure 
including a voxel of uiterest and its selected neighbors; 

Figure 1 8 A is an exemplary image slice from a CT scan of a human 
abdominal region, primarily illustrating a region including the lungs; 
15 Figure 1 8B is a pictorial diagram illustrating the identification of the 

hmg region in the image shce of Figure 1 8 A; 

Figure 1 8C is a pictorial diagram illustrating the removal of the lung 
volume identified in Figure 1 8B; 

« 

Figure 19A is a exemplary image slice form a CT scan of a human 
10 abdominal region, primarily illustrating a region including a portion of the colon and 
bone; 
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Figure 19B is a pictorial diagram illustrating the identificatioii of the 
colon and bone region from the image slice of Figure 1 9 A; 

Figure 19C is a pictorial diagram illustrating the image scan of figure 
1 9a with the regions of bone removed; and 
5 Figure 20 is a flowchart illustrating a method for applying texture to 

monochrome image data 

Figure 21 is a flowchart illustrating a method for volume rendering 
employing a fast perspective ray casting technique; 

Figure 22 is a flowchart illustrating a method for determining the 
1 0 central flight path through a colon lumen employing a volume shrinking technique. 

Figure 23 is a flowchart flirther illustrating a volume shrinking 
technique for use in the method illustrated in Figure 22. 

Figure 24 is a three dimensional pictorial representation of a 
segmented colon lumen with a central fly-path generated therein. 
1 S Figure 25 is a flow chart illustrating a method of generatmg a central 

flight path tbrough a colon lumen employing a segmentation technique. 

Figure 26 is a block diagram of a system embodiment based on a 
personal computer bus architecture. 

Figure 27 is a flow chart illustrating a method of performing volume 
20 imaging using the system of Figure 26. 

Figure 28 is a flow chart illustrating a multi-scan method for 
performing virtual examination of an object, such as a bladder (virtual cystoscopy). 

Figure 29 is a pictorial representation of a display window suitable for 
presenting imaging results from the virtual cystoscopy method of Figure 28 and 
25 providing illustrative outside views of a bladder structure. • 

Figure 30 is a pictorial representation of a display window suitable for 
presenting imaging results from the virtual cystoscopy method of Figure 28 and 
providing illustrative interior views of a bladder structure. 

Figure 31 is a flow chart of a method of performing virtual 
30 examination of an object, such as the larynx, using multiresolution viewing. 
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Figure 32 is a flow chart of a method for performing virtual 

angiography. 

Figures 33 A-C are pictorial views of a portion of the aorta illustrating 
the presence of an abdominal aortic aneurysm. 
. 5 Figure 34 is a flow chart illustrating a method for generating a skeleton 

structure of an object. 

Figure 35 is a schematic diagram of a 26-connected, Euclidean 
weighted, cubic distance plate. 

Figure 36 is a diagram illustrating pseudo-code of a process for 
10 generating a distance map for use in the method of Figure 34. 

DETAILED DESCRIPnQN QF PREFERRED EMBODIMENTS 

While the methods and systems described in this application can be 
applied to any object to be examined, the preferred embodiment which wUl be 
described is the examination of an organ in the human body, specifically the colon. 
15 The colon is long and twisted which makes it especially suited for a virtual 

examination saving the patient both money and the discomfort and danger of a 
physical probe. Other examples of organs which can be examined, without limitation, 
include the lungs, stomach and portions of the gastro-intestinal system, the heart and 
bloodvessels. 

20 Fig. 1 illustrates the steps necessary to perform a virtual colonoscopy 

using volume visualization techniques. Step 101 prepares the colon to be scanned in 
order to be viewed for examination if required by either the doctor or the particular 
scanning instrument This preparation could include cleansing the colon with a 
"cocktail" or liquid which enters the colon after being orally ingested and passed 

25 through the stomach. The cocktail forces the patient to expel waste material that is 
present in the coloiL One example of a substance used is Golytely. Additionally, in 
the case of the colon, air or CO2 can be forced into the colon in order to expand it to 
make the colon easier to scan and examine. This is accomplished with a small tube 
placed in the rectum with approximately 1 ,000 cc of air pumped into the colon to 

30 distend the colon. Depending upon the type of scanner used, it may be necessary for 
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the patient to drink a contrast substance such as barium to coat any unexpunged stool 
in order to distmguish the waste in the colon from the colon walls themselves. 
Altematively, the method for virtually examining the colon can remove the virtual 
waste prior to or during the virtual examination as ejqplained later in this specification. 
5 Step 101 does not need to be perfonned in all examinations as indicated by the dashed 
hneinFig. 1. 

Step 103 scans the organ which is to be examined. The scanner can be 
an apparatus well known in the art, such as a spiral CT-scanner for scanning a colon 
or a Zenita MRI machine for scanning a lung labeled for example with xenon gas. 

10 The scanner must be able to take multiple images from different positions around the 
body during suspended respiration, in order to produce the data necessaiy for the 
volume visualization. An example of a single CT-image would use an X-ray beam of 
5mm width, 1 : 1 to 2: 1 pitch, with a 40cm field-of-view being performed from the top 
of the splenic flexure of the colon to the rectum. 

1 5 Discrete data representations of said obj ect can be produced by other 

methods besides scanning. Voxel data representing an object can be derived from a 
geometric model by techniques described in U.S. Pat. No. 5,038,302 entitled "Method 
of Converting Continuous Three-Dimensional Geometrical Representations into 
Discrete Three-Dimensional Voxel-Based Representations Within a Three- 

20 Dimensional Voxel-Based System" by Kaufinan, issued Aug. 8, 1991, filed My 26, 
1988, which is hereby incorporated by reference. Additionally, data can be produced 
by a computer model of an image which can be converted to three-dimension voxels 
■ and explored in accordance with this invention. One example of this type of data is a 
computer simulation of the turbulence surrounding a space shuttle craft. 

25 Step 104 converts the scanned images into three-dimensional volume 

elements (V oxels). In the preferred embodiment for examining a colon, the scan data 
is reformatted into 5mm thick slices at increments of 1mm or 2.5mm and 
reconstructed in 1mm slices, with each sUce represented as a matrix of 512 by 
512 pixels. By doing this, voxels of approximately 1 cubic mm are created. Thus a 

30 large number of 2D slices are generated depending upon the length of the scan. The 
set of 2D slices is then reconstructed to 3D voxels. The conversion process of 2D 
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images from the seamier into 3D voxels can either be perfonned by the scamiing 
machine itself or by a s^arate machine such as a computer with techniques which are 
well known in the art (for example, see U.S. Pat. No. 4,985,856 entitled "Method and 
Apparatus for Storing, Accessing, and Processing Voxel-based Data" by Kaufm a n et 
5 al.; issued Jan. 15, 1 991, filed Nov. 1 1, 1988; which is hereby incorporated by 
reference). 

Step 105 allows the operator to define the portion of the selected organ 
to be examined. A physician may be interested in a particular section of the colon 
likely to develop polyps. The physician can view a two dimensional slice overview 

10 map to indicate the section to be examined. A starting point and finishing point of a 
path to be viewed can be indicated by the physician/operator. A conventional 
computer and computer interface (e.g., keyboard, mouse or spaceball) can be used to 
designate the portion of the colon which is to be inspected. A grid system with 
coordinates can be used for keyboard entry or the physician/operator can "click" on 

15 the desired points. The entire image of the colon can also be viewed if desired. 

Step 107 performs the plaimed or guided navigation operation of the 
virtual organ being examined. Performing a guided navigation operation is defined as 
navigating through an environment along a predefined or automatically predetermined 
flight path which can be manually adjusted by an operator at any time. After the scan 

20 data has been converted to 3D voxels, the inside of the organ must be traversed from 
the selected start to the selected finishing point The virtual examinations is modeled 
on having a tiny camem traveling through the virtual space with a lens pointing 
towards the finishing point The guided navigation technique provides a level of 
interaction with the camera, so that the camera can navigate through a virtual 

25 environment automatically in the case of no operator interaction, and at the same time, 
allow the operator to manipulate the camera when necessary. The preferred 
embodiment of achieving guided navigation is to use a physically based camera model 
which employs potential fields to control the movement of the camera and which are 
described in detail in Figs. 2 and 3. 

30 Step 109, which can be perfonned concurrently with step 107, displays 

the inside of the organ from the viewpoint of the camera model along the selected 
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. pathway of the guided navigation operation. Three-dimensional displays can be 
generated using techniques well known in the art such as the marching cubes 
technique. However, in order to produce a real time display of the colon, a technique 
is required which reduces the vast number of computations of data necessaiy for the 

5 display of the virtual organ. Fig. 9 describe this display step in more detail. 

The method described in Figure 1 can also be applied to scanning 
multiple organs in a body at the same time. For example, a patient may be examined 
for cancerous growths in both the colon and lungs. The method of Figure 1 would be 
modified to scan all the areas of interest in step 1 03 and to select the current organ to 

10 be examined in step 1 05. For example, the physician/operator may initially select the 
colon to virtually e)q)lore and later explore the linng. Alternatively, two different 
doctors with different specialties may virtually explore different scanned organs 
relating to their respective specialties. Following step 109, the next organ to be 
examined is selected and its portion will be defined and explored. This continues 

1 5 until all organs which need examination have been processed. 

The steps described in conjunction with Figure 1 can also be applied to 
the exploration of any object which can be represented by volume elements. For 
example, an architectural structure or inanimate object can be represented and 
explored in the same manner. 

20 Figure 2 depicts a "submarine" camem control model which performs 

the guided navigation technique in step 107. When there is no operator control during 
guided navigation, the default navigation is similar to that of planned navigation 
which automatically directs the camera along a flight path firom one selected end of 
the colon to another. During the plaimed navigation phase, the camera stays at the 

25 center of the colon for obtaining better views of the colonic surface. When an 
interesting region is encountered, the operator of the virtual camera using guided 
navigation can interactively bring the camera close to a specific region and direct the 
motion and angle of the camera to study the interesting area in detail, without 
unwillingly colliding with the walls of the colon. The operator can control the camera 

30 with a standard interface device such as a keyboard, mouse or non-standard device 
such as a spaceball. In order to fiiUy operate a camera in a virtual environment, six 
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degrees of freedom for the camera is required. The camera must be able to move in 
the hojcizontal, vertical, and Z direction (axes 217), as well as being able to rotate in 
another three degrees of freedom (axes 219) to allow the camera to move and scan all 
sides and angles of a virtual environment The camera model for guided navigation 
includes an inextensible, weightless rod 201 coxmecting two particles 203 and Xj 
203, both particles being subjected to a potential freld 215. The potential freld is 
defined to be highest at the walls of the organ in order to push the camera away from 
the walls. 

The positions of the particles are given by Xj and Xj, and they are 
assumed to have the same mass m. A camera is attached at the head of the submarine 
Xi 203, whose viewing direction coiacides with XjXp The submarine can perform . 
translation and rotation around the center of mass x of the model as the two particles 
are affected by the forces from the potential field V(x) which is defined below, any 
friction forces, and any simulated external force. The relations between Xj, Xj, and x 
are as follows: 

X = (Xjy,z)y 

r = (rsin0cos^,rsinOsin<p,rcos0), 
Xj = x + r, 

X2^x-r, (1) 

where r, 9 and cp are the polar coordinates of the vector xxj. 

The kinetic energy of the model, T, is defined as the summation of the kinetic 

energies of the movements of Xi and X2: 

• 2 I '2 

= mx +mp'* 

= m{x^ + y^ + z^) + mr^{6^ + ^hm^). (2) 

Then, the equations for the motion of the submarine model are 
obtained by using LaGrange's equation: 
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(3) 



where the qjS are the generalized coordinates of tibie model and can be considered as 
the variables of time / as: 



with \|/ denoting the roll angle of our camera system, which will be explained later. 
The FjS axe called the generalized forces. The control of the submarine is performed 
5 by applying a simulated external force to x,. 



and it is assumed that both and Xj ar6 affected by the forces from the potential field 
and the frictions which act in the opposite direction of each particle's velocity. Conse- 
quently, the generalized forces are formulated as follows: 



1 0 where k denotes the friction coefficient of the system. The external force F^^ is 

applied by the operator by simply clicking the mouse button in the deshed direction 
207 in the generated image, as shown in Figure 2. This camera model would then be 
moved in that direction. This allows the operator to control at least five degrees of 
freedom of the camera with only a smgle click of the mouse button. From Equations 

1 5 (2), (3) and (5), it can be derived that the accelerations of the five parameters of our 
submarine model as: 



(4) 




Fj. = -mVy(x2) - fcx2, 



est, 



(5) 
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y = 



2^ dx dx ^ m 2m' 

2^ dy ^ By ^ m 2m' 
5 , ^( dVjxr) , flV(x2) , kz . F, 
^ ^ 2^ dz dz ^ m"*'2m' 

[cos «{co.^^ - ^) + sin - ^)} 

^ dz dz 

Jfc • 1 

0 + - (F,cCO3Bcos^ + FjfCos0sm<f>— FgSmB)^ 

- Tlx 277ir 



--{- sm <K—^ ^) + cos<Pi-g^ Q^)} 

^siatf + :r^(--fa. sin ^ + fy cos ^)], (6) 



where x and x denote the first and the second derivative of x, respectively, and 



[ ^YS^ , dV(xl ^ dV(x) \ jjgjjotgg tijg gradient of the potential at a point x. 
^ ax dy dz ) 



Theterms ^^sinecosQ of 6 and -i£l£2£i of ^" are called the 

sine 



centrifugal force and the Coriolis force^ respectively, and they are concerned with the 
5 exchange of angular velocities of the submarine. Since the model does not have the 
moment of inertia defined for the rod of the submarine, these terms tend to cause an 
overflow of the numeric calculation of q>. Fortunately, these terms become significant 
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only when the angular velocities of the submarine model are significant, which 
essentially means that the camera moves too fast Since it is meaningless to allow the 
camera to move so fast because the organ could not be properly viewed, these terms 
are minimized in our implementation to avoid the overflow problem. 
) From the first three formulas of Equation (6), it is known that the 

submarine cannot be propelled by the external force against the potential field if the 
following condition is satisfied: 



|vnxO + vnxOI> 



m 



Since the velocity of the submarine and the external force have upper limits in our 

implementation, by assigning sufBciently high potential values at the boundary of the 
1 0 objects, it can be guaranteed that the submarine never bumps against the objects or 

walls in the environment 

As mentioned previously, the roll angle \|/ of the camera system needs 

to be considered. One possible option' allows the operator full control of the angle \|/. 

However, although the operator can rotate the camera fireely around the rod of the 
15 model, he or she can easily become disoriented. The preferred technique assumes that 

the upper direction of the camera is connected to a pendulimi with mass m2 U 

which rotates fireely around the rod of the submarine, as shown in Figure 3. The 

direction of the pendulum, r2, is expressed as: 



rj — r2(cos 0 cos ip sm ^ + sin ^ cos ^, cos sin 0 sin -0 — cos <f> cos — sin 6 sin -^j. 



although it is possible to calculate the accurate movement of this pendulum along 
20 with the movement of the submarine, it makes the system equations too complicated. 
Therefore, it is assumed that all the generalized coordinates except the roll angle \|/ are 
constants, and thus define the independent kinetic energy for the pendulum system as: 
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m2 .2 m2r2^ ]2 



This simplifies the model for the roll angle. Since it is assumed in this model that the 
gravitational force 

Fy = m2g = (Wl25a:,m2<7y,m2^x) 



acts at the mass point m^i the acceleration of \|/ can be derived using LaGrange's 
equation as: 

(cos d sin ^ cos -0 + cos ^ sin ^) 

sin^cos^)} (7) 

7112 



5 From Equations (6) and (7), the generalized 



coordinates q(t) and their derivatives q(t) are calculated asymptotically by using 
Taylor series as: 

q(i + M = q(f) + /iq(t) + 0(fe'), 



10 



to &eely move the submarine. To smooth the submarine's motion, the time step h is 
selected as an- equilibrimn value betwera being as small as possible to smooth the 
motion but as large as necessary to reduce computation cost. 
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Definition of the Potential Field 
The potential field in the submarine model in Figure 2 defines the 
boundaries (walls or other matter) in the virtual organ by assigning a high potential to 
the boundary in order to ensure that the submarine camera does not collide with the 
5 walls or other boundary. If the camera model is attempted to be moved into a high 
potential area by the operator, the camera model will be restrained firom doing so 
unless the operator wishes to examine the organ behind the boxmdary or inside a 
polyp, for example. In the case of performing a virtual colonoscopy, a potential field 
value is assigned to each piece of volumetric colon data (volume element). When a 

1 0 particular region of interest is designated in step 105 of Fig. 1 with a start and finish 
point, the voxels within the selected area of the seamed colon are identified using 
conventional blocking operations. Subsequently, a potential value is assigned to every 
voxel X of the selected volume based on the following three distance values: the 
distance firom the finishing point dt(x), the distance firom the colon surface ds(x) and 

15 the distance firom the center-line of the colon space dc(x). dt(x) is calculated by usmg 
a conventional growing strategy. The distance firom the colon surface, ds(x), is 
computed using a conventional technique of growing from the surface voxels inwards. 
To determine dc(x), the center-line of the colon from the voxel is first extracted, and 
then dc(x) is computed usiag the conventional growing strategy from the center-line 

20 of the colon. 

To calculate the center-line of the selected colon area defined by the 
user-specified start point and the user-specified finish point, the maximum value of 
ds(x) is located and denoted dmax. Then for each voxel inside the area of interest, a 

« 

cost value of dmax - ds(x) is assigned. Thus the voxels which are close to the colon 
25 surface have high cost values and the voxels close to the center line have relatively 
low cost values. Then, based on the cost assignment, the single-source shortest path 
technique which is well known in the art is applied to ejBBciently compute a minimn Tn 
cost path from the source point to the finish point. This low cost line iadicates the 
center-line or skeleton of the colon section which is desired to be explored. This 
30 technique for determining the center-line is the preferred technique of the invention. 
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To compute the potential value V(x) for a voxel x inside the area of 
interest, the following formula is employed: 



V{x) = C^dt{^r + C2{ 




(8) 



4(x) + d,{x) 



where Cj, C2, \i and v are constants chosen for the task. In order to avoid any 
collision between the virtual camera and the virtual colonic surface, a sufGciently 
5 large potential value is assigned for all points outside the colon. The gradient of the 
potential field will therefore become so significant that the submarine model camera 
will never collide with the colonic wall when being run. 



two side wails 401 and 403 of the colon being shown. Two blocldng walls are 
selected by the operator in order to define the section of the colon which is of interest 
to examine. Nothing can be viewed beyond the blocking walls. This helps reduce the 
number of computations when displaying the virtual representation. The blocking 
15 walls together with side walls identify a contained volumetric shape of the colon 
which is to be explored. 



examination, the start volume element 501 and the finish volume element 503. The 
start and finish points are selected by the operator in step 105 of Fig. 1 . The voxels 
20 between the start and finish points and the colon sides are identified and marked, as 
indicated by the area designated with "x"s in Fig. 6. The voxels are three-dimensional 
representations of the picture element. 



voxels in Fig. 6. The outermost layer of all the voxels (closest to the colon walls) is 
25 peeled off step-by-step, until there is only one inner layer of voxels remaining. Stated 
differently, each voxel furthest away from a center point is removed if the removal 
does not lead to a disconnection of the path between the start voxel and the finish 



10 



Another technique to determine the center-line of the path in the colon 
is called the "peeHayer" technique and is shown in Figure 4 through Figure 8. 

Figure 4 shows a 2D cross-section of the volxmietric colon, with the 



Figure 5 shows two end points of the flight path of the virtual 



The peel-layer technique is then applied to the identified and marked 
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voxel. Figure 7 shows the intemediate result after a number of iterations of peehng 
the voxels in the virtual colon are complete. The voxels closest to the walls of the 
colon have been removed. Fig. 8 shows the final flight path for the camera model 
down the center of the colon after all the peeling iterations are complete. This 
5 produces essentially a skeleton at the center of the colon and becomes the desired 
flight path for the camera model. 

Ir Buffer assisted visibility 
Figure 9 describes a real time visibility technique to display of virtual 
images seen by the camera model in the virtual three-dimensional volume 

1 0 representation of an organ. Figure 9 shows a display technique using a modified Z 
buffer which corresponds to step 109 in Fig. L The number of voxels which could be 
possibly viewed firom the camera model is extremely large. Unless the total number 
of elements (or polygons) which must be computed and visualized is reduced firom an 
entire set of voxels in the scanned environment, the overall number of computations 

1 5 will make the visualization display process exceedingly slow for a large internal area. 
However, in the present invention only those images which are visible on the colon 
surface need to be computed for display. The scanned environment can be subdivided 
into smaller sections, or cells. The Z buffer technique then renders only a portion of 
the cells which are visible firom the camera. The Z buffer technique is also used for 

20 three-dimensional voxel representations. The use of a modified Z buffer reduces the 
number of visible voxels to be computed and allows for the real time examination of 
the virtual colon by a physician or medical technician. 

The area of interest from which the center-line has been calculated in 
step 107 is subdivided into cells before the display technique is apphed. Cells are 

25 collective groups of voxels which become a visibility unit. The voxels m each cell 
will be displayed as a group. Each cell contains a number of portals through which 
the other cells can be viewed. The colon is subdivided by beginning at the selected 
start point and moving along the center-line 1 00 1 towards the finish point. The colon 
is then partitioned into cells (for example, cells 1003, 1005 and 1007 in Fig. 10) when 

30 a predefined threshold distance along the center-path is reached. The threshold 
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distance is based i^on the specifications of the platform upon which the visualization 
technique is performed and its capabilities of storage and processing. The cell size is 
directly related to the nximber of voxels which can be stored and processed by the 
platform. One example of a threshold distance is 5cm, although the distance can 
5 greatly vary. Each cell has two cross-sections as portals for viewing outside of the 

■ 

cell as shown in Fig. 10. 

Step 901 in Fig, 9 identifies the cell within the selected organ which 
currently contains the camera. The current cell will be displayed as well as all other 
cells which are visible given flie orientation of the camera. Step 903 builds a stab tree 

10 (tree diagram) of hierarchical data of potentially visible cells from the camera 

(through defined portals), as will be described in fiirther detail hereinbelow. The stab 
tree contains a node for every cell which may be visible to the camera. Some of the 
cells may be transparent without any blocking bodies present so that more than one 
cell will be visible in a smgle direction. Step 905 stores a subset of the voxels from a 

1 5 cell which include the intersection of adjoining cell edges and stores them at the 
outside edge of the stab tree in order to more efficiently determine which cells are 
visible. 

Step 907 checks if any loop nodes are present in the stab tree. A loop 
node occurs when two or more edges of a single cell both border on the same nearby 
20 cell. This may occur when a single cell is surrounded by another cell. If a loop node 
is identified in the stab tree, the method continues with step 909. If there is no loop 
node, the process goes to step 911. 

Step 909 collapses the two cells making up the loop node into one 
large node. The stab tree is then corrected accordingly. This eliminates the problem 
25 of viewing the same cell twice because of a loop node. The step is performed on all 
identified loop nodes. The process then continues with step 911. 

Step 91 1 then mitiates the Z-buffer with the largest Z value. The Z 
value defines the distance away from the camera along the skeleton path. The tree is 
then traversed to first check the intersection values at each node. If a node 
30 intersection is covered, meaning that the current portal sequence is occluded (which is 
detenmned by the Z buffer test), then the traversal of the cunent branch in the tree is 
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Stopped. Step 913 traverses each of the branches to check if the nodes are covered 
and displays them if they are not. 

Step 915 then constructs the image to be displayed on the operator's 
screen firom the volume elements within the visible cells identified in step 913 using 
5 one of a variety of techniques known in the art, such as volume rendering by 
compositing. The only cells shown are those which are identified as potentially 
visible. This technique limits the number of cells which reqidres calcxilations in order 
to achieve a real time display and conespondingly increases the speed of the display 
for better performance. This technique is an improvement over prior techniques 
which calculate all the possible visible data points whether or not they are actually 
viewed. 

Figure 11 A is a two dimensional pictorial representation of an organ 
which is being explored by guided navigation and needs to be displayed to an 
operator. Organ 1101 shows two side walls 1 102 and an object 1 105 in the center of 
the pathway. The organ has been divided into four cells A 1 151, B 1 153, C 1 155 and 
D 1157. The camera 1 103 is facing towards cell D 1157 and has a field of vision 
defined by vision vectors 1 107, 1 108 which can identify a cone-shaped field. The 
cells which can be potentially viewed are cells B 1 1 53, C 1 1 55 and D 1 1 57. Cell C 
1 155 is completely surrounded by Cell B and thus constitutes a node loop. 

Fig. 1 IB is a representation of a stab tree built firom the cells in Fig. 
1 1 A. Node A 1 1 09 which contains the camera is at the root of the tree. A sight line 
or sight cone, which is a visible path without being blocked, is drawn to node B 1 1 10. 
Node B has direct visible sight lines to both node C 1 1 12 and node D 1 1 14 and which 
is shown by the connecting arrows. The sight line of node C 1 1 12 in the direction of 
the viewing camera combines with node B 1 1 10. Node C 1 1 12 and node B 1 1 10 will 
thus be collapsed hito one large node B' 1122 as shown in Fig. 1 IC. 

Fig. lie shows node A 1 109 containing the camera adjacent to node B' 
1 122 (containing both nodes B and node C) and node D 1 1 14. The nodes A, B' and D 
will be displayed at least partially to the operator. 

Figs 12A - 12E illustrate the use of the modified Z buffer with cells 
that contain objects which obstruct the views. An object could be some waste 
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material in a portion of the virtual colon. Fig. 12A shows a virtual space with 
10 potential cells: A 1251, B 1253, C 1255, D 1257, E 1259, F 1261, G 1263, H 1265, 
1 1267 and J 1269. Some of the cells contain objects. If the camera 1201 is 
positioned in cell 1 1267 and is facing toward cell F 1261 as indicated by the vision 

5 vectors 1203, then a stab tree is generated in accordance with the technique illustrated 
by the flow diagram in Fig. 9. Fig. 12B shows the stab tree generated with the - 
intersection nodes showing for the virtual representation as shown in Fig. 12 A. Fig. 
12B shows cell 1 1267 as the root node of the tree because it contains the camera 
1201. Node 1 121 1 is pointing to node F 1213 (as mdicated with an arrow), because 

10 cell F is directly connected to the sight line of the camera. Node F 1213 is pointing to 
both node B 1215 and node E 1219. Node B 1215 is pointing to node A 1217. Node 
C 1202 is completely blocked from the line of sight by camera 1201 so is not included 
in the stab tree. 

Fig. 12C shows the stab tree after node 1 121 1 is rendered on the 
1 5 display for the operator. Node 1 1 2 1 1 is then removed from the stab tree because it 
has already been displayed and node F 1213 becomes the root Fig. 12D shows that 
node F 1213 is now rendered to join node 11211. The next nodes in the tree 
connected by arrows are then checked to see if they are aheady covered (already 
processed). In this example, all of the intersected nodes from the camera positioned 
20 in cell 1 1 267 has been covered so that node B 5 1 5 (and therefore dependent node A) 
do not need to be rendered on the display. 

Fig. 12E shows node E 515 being checked to determine if its 
intersection has been covered. Since it has, the only rendered nodes in this example 
of Figure 12A-12E are nodes I and F while nodes A, B and E are not visible and do 
25 not need to have their cells prepared to be displayed. 

The modified Z buffer technique described in Figure 9 allows for fewer 
computations and can be applied to an object which has been represented by voxels or 
other data elements, such as polygons. 

Figure 13 shows a two dimensional virtual view of a colon with a large 
30 polyp present along one of its walls. Figure 13 shows a selected section of a patient's 
colon which is to be examined ftuther. The view shows two colon walls 1301 and 
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1303 with the growth indicated as 1305. Layers 1307, 1309, and 1311 show inner 
layers of the growflL It is desirable for a physician to be able to peel the layers of the 
polyp or tumor away to look inside of the mass for any cancerous or other harmful 
material. This process would in effect perform a virtual biopsy of the mass without 

5 actually cutting into the mass. Once the colon is represented virtually by voxels, the 
process of peeling away layers of an object is easily performed in a similar manner as 
described in conjunction with Figs. 4 tlirough 8. The mass can also be sliced so that a 
particular cross-section can be examined. In Fig. 13, a planar cut 13 13 can be made 
so that a particular portion of the growth can be examined. Additionally, a user- 

0 defined slice 1319 can be made in any manner in the growth. The voxels 1319 can 
either be peeled away or modified as explained below. 

A transfer fiinction can be performed to each voxel in the area of 
interest which can make the object transparent, semi-transparent or opaque by altering 
' coefficients representing the translucently for each voxel. An opacity coefficient is 

5 assigned to each voxel based on its density. A mapping function then transforms the 
density value to a coefficient representing its translucency. A high density scanned 
voxel will indicate either a wall or other dense matter besides simply open space. An 
operator or program routine could then change the opacity coefficient of a voxel or 
group of voxels to make them appear transparent or semi-transparent to the submarine 

0 camera model. For example, an operator may view a tumor within or outside of an 
entire growth. Or a transparent voxel will be made to appear as if it is not present for 
the display step of Figure 9. A composite of a section of the object can be created 
using a weighted average of the opacity coefficients of the voxels in that section. 

If a physician desires to view the various layers of a polyp to look for a 

15 cancerous areas, this can be performed by removing the outer layer of polyp 1305 
yielding a first layer 1307. Additionally, the first irmer layer 1307 can be stripped 
back to view second inner layer 1309. The second inner layer can be stripped back to 
view third inner layer 131 1, etc. The physician could also slice the polyp 1305 and 
view only those voxels within a desired section. The slicing area can be completely 

10 user-defined. 
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Adding an opacity coefficient can also be used in other ways to aid in 
the exploration of a virtual system. If waste material is present and has a density as 
other properties within a certain known range, the waste can be made transparent to 
the virtual camera by changing its opacity coefBcient during the examination. This 

5 will allow the patient to avoid ingesting a bowel cleansing agent before the procedure 
and make the examination faster and easier. Other objects can be similarly made to 
disappear dq)ending upon the actual application. Additionally, some objects like 
polyps could be enhanced electronically by a contrast agent followed by a use of an 
' appropriate transfer function. 

0 Figure 14 shows a system for performing the virtual examination of an 

object such as a human organ using the techniques described in this specification. 
Patient 1401 lies down on a platform 1402 while scanning device 1405 scans the area 
that contains the organ or organs which are to be examined. The scanning device 
1405 contains a scanning portion 1403 which actually takes images of the patient and 

5 an electrojoics portion 1406. Electronics portion 1406 comprises an interface 1407, a 
central processing unit 1409, a memory 141 1 for temporarily storing the scanning 
data, and a second interface 1413 for sending data to the virtual navigation platform. 
Interface 1407 and 1413 could be included in a single interface component or could be 
the same component. The components in portion 1406 are connected together with 

0 conventional connectors. 

In system 1400, the data provided from the scanning portion of device 
1403 is transferred to portion 1405 for processing and is stored in memory 1411. 
Central processing unit 1409 converts the scanned 2D data to 3D voxel data and 
stores the results in another portion of memory 141 1 . Alternatively, the converted 

5 data could be directly sent to interface unit 141 3 to be transferred to the virtual 

navigation terminal 1416. The conversion of the 2D data could also take place at the 
virtual navigation terminal 1416 after being transmitted from interface 1413. In the 
preferred embodimeat, the converted data is traosmitted over carrier 1414 to the 
virtual navigation terminal 1416 in order for an operator to perform the virtual 

0 examination. The data could also be transported in other conventional ways such as 
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Storing the data on a storage medium and physically transporting it to terminal 1416 or 
by using satellite transmissions. 

The scanned data may not be converted to its 3D representation imtil 
the visualization rendering engine requires it to be in 3D form; This saves 
5 computational steps and memory storage space. 

Virtual navigation terminal 1416 includes a screen for viewing the 
virtual organ or other scaimed image, an electronics portion 1415 and interface control 
1419 such as a keyboard, mouse or spaceball. Electronics portion 1415 comprises a 
interface port 1421, a central processing unit 1423, other components 1427 necessary 

10 to run the terminal and a memory 1425. The components in terminal 1416 are 
connected together with conventional connectors. The converted voxel data is 
received in interface port 1421 and stored in memory 1425. The central processor 
unit 1423 then assembles the 3D voxels into a virtual representation and runs the 
submarine camera model as described in Figures 2 and 3 to perform the virtual 

15 examination. As the subrnarine camera travels through the virtual organ, the visibility 
technique as described in Figure 9 is used to compute only those areas which are 
visible from the virtual camera and displays them on screen 1417. A graphics 
accelerator can also be used in generating the representations. The operator can use 
interface device 1419 to indicate which portion of the scaimed body is desired to b,e 

20 explored. The interface device 1419 can further be used to control and move the 
submarine camera as desired as discussed in Figure 2 and its accompanying 
description. Terminal portion 1415 can be the Cube-4 dedicated system box, 
generally available from the Department of Computer Science at the State University 
of New York at Stony Brook. 

25 Scanning device 1405 and terminal 1416, or parts thereof, can be part 

of the same unit A single platform would be used to receive the scan image data, 
connect it to 3D voxels if necessary and perform the guided navigation. 

An important feature in system 1400 is that the virtual organ can be 
examined at a later time without the presence of the patient. Additionally, the virtual 

30 examination could take place while the patient is being scanned. The scan data can 
also be sent to multiple terminals which would allow more than one doctor to view 
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the inside of the organ simultaneously. Thus a doctor in New York could be looking 
at the same portion of a patient's organ at the same time with a doctor in California 
while discussing the case. Alternatively, the data can be viewed at different times. 
Two or more doctors could perform their own examination of the same data in a 
difficult case. Multiple virtual navigation terminals could be used to view the same 
scan data. By reproducing the organ as a virtual organ with a discrete set of data, 
there are a multitude of benefits in areas such as accuracy, cost and possible data 
manipiQations. 

The above described techniques can be further enhanced in virtual 
colonoscopy applications through the use of an improved electronic colon cleansing 
technique which employs modified bowel preparation operations followed by image 
segmentation operations, such that fluid and stool remaining in the colon during a 
computed tomographic (CT) or magnetic resonance imaging (MRI) scan can be 
detected and removed fi:om the virtual colonoscopy images. Through the use of such 
techniques, conventional physical washing of the colon, and its associated 
inconvenience and discomfort, is minimized or completely avoided. 

Referring to Figure 15, the first step ia electronic colon cleansing is 
bowel preparation (step 1510), which takes place prior to conducting the CT or 
magnetic resonance imaging (MRI) scan and is intended to create a condition where 
residual stool and fluid remaining in the colon present significanfly different image 
properties from that of the gas-filled colon interior and colon wall. An exemplary 
bowel preparation operation includes ingesting three 250 cc doses of Barium Sulfate 
suspension of 2,1 % WA^, such as manufactured by E-Z-EM, Inc.,of Westbury, New 
York, during the day prior the CT or MRI scan. The three doses should be spread out 
over the course of the day and can be ingested along with three meals, respectively. 
The Barium Sulfate serves to enhance the images of any stool which remains in the 
colon. In addition to the intake of Barium Sulfate, fluid intake is preferably mcreased 
dxiring the day prior to the CT or MRI scan. Cranberry juice is known to provide 
increased bowel fluids and is preferred, although water can also be ingested. In both 
30 the evening prior to the CT scan and the morning of the CT scan, 60 .ml of a 

Diatrizoate Meglumine and Diaztrizoate Sodium Solution, which is commercially 
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available as MD-Gastroview, manufactured by Mallinckrodt, lac. of St Louis, 
Missouri, can be consumed to enhance image properties of the colonic fluid. Sodium 
phosphate can also be added to the solution to liquidize the stool in the colon, which 
provides for more uniform enhaacement of the colonic fluid and residual stool. 
5 The above described exemplary preliminary bowel preparation 

operation can obviate the need for conventional colonic washing protocols, which can 
call for the ingestion of a gallon of Golytely solution prior to a CT scan. 

Just prior to conducting the CT scan, an intravenous injection of 1 ml 
of Glucagon, manufactured by Ely Lily and Company, of Indianapohs, Indiana can be 

1 0 administered to minimize colon collapse. Then, the colon can be inflated using 
approximately lOOOcc of compressed gas, such as COj, or room air, which can be 
introduced through a rectum tube. At this point, a conventional CT scan is performed 
to acquire data from the region of the colon (step 1 520). For example, data can be 
acquired using a GE/CTI spiral mode scanner operating in a helical mode of Smm, 

15 1.5-2.0:1 pitch, reconstructed in 1mm slices, where the pitch is adjusted based upon 
the patient's height in a known mamer. A routine imaging protocol of 120 kVp and 
200-280 ma can be utilized for this operation. The data can be acquired and 
reconstructed as 1mm thick shce images having an array size of 512x512 pixels in the 
field of view, which varies from 34 to 40 cm depending on the patient's size, the 

20 number of such slices generally varies under these conditions from 300 to 450, 
depending on the patient's height. The image data set is converted to volume 
elements or voxels (step 1530). 

Image segmentation can be performed in a number of ways. In one 
present method of image segmentation, a local neighbor technique is used to classify 

25 voxels of the image data in accordance with similar intensity values. In this method, 
each voxel of an acquired image is evaluated with respect to a group of neighbor 
voxels. The voxel of interest is referred to as the central voxel and has an associated 
intensity value. A classification indicator for each voxel is established by comparing 
the value of the central voxel to each of its neighbors. If the neighbor has the same 

30 value as the central voxel, the value of the classification indicator is incremented. 
However, if the neighbor has a different value from the central voxel, the 
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classification indicator for the central voxel is decremented. The central voxel is then 
classified to that category which has the maximum indicator value, which indicates 
the most uniform neighborhood among the local neighbors. Each classification is 
indicative of a particular intensity range, which in turn is representative of one or 
more material types being imaged. The method can be further enhanced by 
employing a mixture probability function to the similarity classifications derived. 

An alternate process of image segmentation is performed as two major 
operations: low level processing and high level feature extraction. During low level 
processing, regions outside the body contour are eliminated from fiuther processing 
and voxels within the body contour are roughly categorized in accordance witii well 
defined classes of intensity characteristics. For example, a CT scan of the abdominal 
region generates a data set which tends to exhibit a well defined intensity distribution. 
The graph of Figure 16 illustrates such an intensity distribution as an exemplary 
histogram having four, well defined peaks, 1602, 1604, 1606, 1608, which can be 
classified according to intensity thresholds. 

The voxels of the abdominal CT data set are roughly classified as four 
clusters by intensity thresholds (step 1540). For example, Cluster 1 can include 
voxels whose intensities are below 140. This cluster generally corresponds to the 
lowest density regions within the interior of the gas filled coIoil Cluster 2 can include 
voxels which have intensity values in excess of 2200. These intensity values 
correspond to the enhanced stool and fluid within the colon as well as bone. Cluster 
3 can include voxels with intensities in the range of about 900 to about 1080. This 
intensity range generally represents soft tissues, such as fat and muscle, which are 
unlikely to be associated with the colon. The remaining voxels can then be grouped 
together as cluster 4, which are likely to be associated with the colon wall (including 
mucosa and partial volume mixtures around the colon wall) as well as lung tissue and 
soft bones. 

Clusters 1 and 3 are not particularly valuable in identifying the colon 
wall and, therefore are not subject to substantial processing during image 
segmentation procedures for virtual colonoscopy. The voxels associated with cluster 
2 are important for segregating stool and fluid firom the colon wall and are processed 
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further during the high-level feature extraction operations. Low level processing is 
concentrated on the fourth cluster, which has the highest likelihood of corresponding 
to colon tissue (step 1550). 

For each voxel in the fourth cluster, an intensity vector is generated 

5 using itself and its neighbors. The intensity vector provides an indication of the 
change in intensity in the neighborhood proximate a given voxel. The number of 
neighbor voxels which are used to establish the intensity vector is not critical, but 
involves a tradeoff between processing overhead and accuracy. For example, a simple 
voxel intensity vector can be established with seven (7) voxels, which includes the 

0 voxel of interest, its front and back neighbors, its left and right neighbors and its top 
and bottom neighbors, all surrounding the voxel of interest on three mutually 
perpendicular axes. Figure 17 is a perspective view illustrating an exemplary intensity 
vector in the form of a 25 voxel intensity vector model, which includes the selected 
voxel 1702 as well as its first, second and third order neighbors. The selected voxel 

5 1702 is the central point of this model and is referred to as the fixed voxel. A planar 
slice of voxels, which includes 12 neighbors on the same plane as the fixed voxel, is 
referred to as the fixed slice 1 704. On adjacent planes to the fixed slice are two 
nearest slices 1706, having five voxels each. Adjacent to the first nearest slices 1706 
are two second nearest slices 1708, each having a single voxel. The collection of 

0 intensity vectors for each voxel hi the fourth cluster is referred to as a local vector 
series. 

Because the data set for an abdominal image generally includes more 
than 300 sHce images, each with a 512 x 5 12 voxel array, and each voxel having an 
associated 25 voxel local vector, it is desirable to perform feature analysis (step 1570) 

5 on the local vector series to reduce the computational burden. One such feature 
analysis is a principal component analysis (PCA), which can be applied to the local 
vector series to determine the dimension of a feature vector series and an orthogonal 
transformation matrix for the voxels of cluster 4- ' 

It has been found that the histogram (Figure 16) of the CT image 

0 intensities tends to be fairly constant from patient to patient for a particular scanner, 
given equivalent preparation and scaiming parameters. Relying on this observation, 
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an orthogonal transfonnation matrix can be established which is a predetermined 
matrix detennined by usmg several sets of training data acquired using the same 
scanner under similar conditions. From this data, a transformation matrix, such as a 
Karlhunen-Lo6ve (Kl-L) transformation, can be generated in a known manner. The 
S transformation matrix is applied to the local vector series to generate feature vector * 
series. Once in the feature-vector space domain, vector quantization techniques can 
be used to classify the feature vector series. 

An analytical, self-adaptive algorithm can be used for the classification 
of the feature vectors. In defining this algorithm, let {XjeR'^ti = 1,2,3,...,N} be the 
10 series of the feature vectors, where N is the number of feature vectors; K denotes the 
maximum number of classes; and T is a threshold which is adaptive to the data set. 
For each class, a representative element is generated by the algorithm. Let a^ be a 
rq)resentative element of class k and n^^ be the number of feature vectors in that class. 
The algorithm can then be outlmed as: 

15 1. Set a^=X^;K^l; 

m 

2. obtain the class number K and class parameters (aj^ 
for(j=l;i<N;i-H-) 

for (j=l; j<K;j++.) 

calculate dj = dist{X^,aj) ; 

20 end for 

index=axc min djj 

i£ {{d^^<T)ior{K=K)) 

update class parameters: 



oc ^ index " "I 7T ^ ^ ^index^ ^ index -^i ) ' 

^index'^index'^^' 
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end if 
else 

generate new class 



^ic + 1 ~ ' 
K = 



end else 



end for 



3. label each feature vector to a class according to the nearest neighbor rule 
for (i = l; i<J^; i + + ) 

. for {i = l; j<K; J+ + ) 
10 calculate = dist {Xj^, a^) ; 

end for 

index ^ arc mm dj. 

label voxel i to class index. 

end for 



15 In this algorithm, dist(x,y) is the Euclidean distance between vector x 

and j; and arc min dj gives the integer / which realizes the minim um value of dy. 

The above described algorithm is dependent only on the parameters T 
and K. However, the value of K, which relates to the number of classes within each 
voxel cluster, is not critical and can be set to a constant value, such as K=l 8. 

20 However, T, which is the vector sjmilarily threshold, greatly influences the 

classification results. If the selected value of T is too large, only a single class will be 
generated. On the other hand, if the value of T is too small, the resulting classes will 
exhibit undesirable redundancy. By setting the value of T to be equal to the maximum 
component variance of the feature vector series, the maximum number of distinct 

25 classes results. 
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As a result of the initial classification process, each voxel within the 
selected cluster is assigned to a class (step 1570). In the exemplary case of virtual 
colonoscopy, there are several classes within cluster 4. Thus, the next task is to 
determine which of the several classes in cluster 4 corresponds to the colon wall. The 

5 first coordinate of the feature vector, which is that coordinate of the ifeature vector . 
exhibiting the highest variance, reflects the information of the average of the 3D local 
voxel intensities. The remaii^g coordinates of the feature vector contaui the 
information of directional intensity change within the local neighbors. Because the 
colon wall voxels for the interior of the colon are generally in close proximity to the 

10 gas voxels of cluster 1 , a threshold interval can be determined by data samples 
selected fi-om typical colon wall intensities of a typical CT data set to roughly 
distinguish colon wall voxel candidates. The particular threshold value is selected for 
each particular imaging protocol and device. This threshold interval can then applied 
to all CT data sets (acquired from the same machine, xxsing the same imaging 

1 5 protocol). If the first coordinate of the representative element is located in the 

threshold interval, the corresponding class is regarded as the colon wall class and all 
voxels ia that class are labeled as colon wall-like voxels. 

Each colon wall-like voxel is a candidate to be a colon wall voxel. 
There are three possible outcomes of not belonging to the colon wall. The first case 

20 relates to voxels which are close to the stool/liquid inside the colon. The second case 

occurs when voxels are in the lung tissue regions. The third case represents mucosa 

* 

voxels. Clearly then, low level classification carries a degree of classification 
uncertainty. The causes of the low-level classification uncertainty vary. For example, 
a partial-volume effect resulting from voxels containing more than one material type 

25 (i.e., fluid and colon wall) leads to the first case of uncertainty. The second and the 
third cases of uncertainty are due to both the partial volume effect as well as the low 
contrast of CT images. To resolve the uncertainty, additional information is needed. 
Thus, a high-level feature extraction procedure is used in the present method to 
further distinguish candidates for the colon wall from other colon wall-like voxels, 

30 based on a,priori anatomical knowledge of the CT images (step 1580). 
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An initial step of the high-level feature extraction procedure can be to 
eliminate the region of lung tissue from the low-level classification results: Figure 
1 8 A is an exemplary slice image clearly illustrating the lung region 1 802. The lung 
region 1802 is identifiable as a generally contiguous three dimensional volume 

5 enclosed by colon wall-like voxels, as illustrated in Figure 1 8B. Given this 

characteristic, the lung region can be identified using a region growing strategy. The 
first step in this technique is to find a seed voxel withm the region of growing. 
Preferably, the operator performing the CT imagmg scan sets the imaging range such 
that the top most slice of the CT scan does not contain any colon voxels. As the 

1 0 interior of lung should be filled with air, the seed is provided by the low-level 

classification simply by selecting an air voxel. Once the lung region outline of Figure 
18B is determined, the limg volume can be removed from the image slice (Figure 
18C). 

A next step hi performing high-level feature extraction can be to 
15 separate the bone voxels from enhanced stool/fluid voxels in cluster 2. The bone 
tissue voxels 1902 are generally relatively far away from the colon wall and resides 
outside the colon volume. To the contrary, the residual stool 1906 and fluid 1904 are 
enclosed inside the colon volume. Combining the a priori proximity information and 
the colon wall mformation obtained from the low-level classification process, a rough 
20 colon wall volume is genemted. Any voxel separated by more than a predetermined 
number (e.g., 3) of voxel units from the colon wall, and outside the colon volume, 
will be labeled as bone and then removed from the image. The remaining voxels in 
cluster 2 can be assumed to represent stool and fluid within the colon volume (see 
Figures 19A-C). 

25 The voxels withhi the colon volume identified as stool 1906 and fluid 

1904 can be removed from the hnage to generate a clean colon lumen and colon wall 
image. In general, there are two kinds of stool/fluid regions. One region type is small 
residual areas of stool 1906 attached to the colon wall. The other region type is large 
volumes of fluid 1904, which collect in basin-like colonic folds (see Figures 19A-C). 

30 The attached residual stool regions 1906 can be identified and removed 

because they are inside the rough colon volume generated during the low-level 
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classification process. The fluid 1906 in the basin-like colon fold usually has a 
horizontal surface 1908 due to the effect of gravity. Above the surface is always a gas 
region, which exhibits a very high contrast to the fluid mtensity. Thus, the surfece 
interface of the fluid regions can be easily marked. 
5 Using a region growing strategy, the contour of the attached stool 

regions 1906 can be outlined, and the part which is away from the colon wall volume 
can be removed. Similarly, the contour of the fluid regions 1904 can also be outlined. 
After eliminating the horizontal surfaces 1908, the colon wall contour is revealed and 

the clean colon wall is obtained. 
10 It is difficult to distinguish the mucosa voxels from the colon wall 

voxels. Even though the above three dimensional processing can remove some 
mucosa voxels, it is difficult to remove all mucosa voxels. In optical colonoscopy, 
physicians directly inspect the colonic mucosa and search for lesions based on the 
color and texture of tihe mucosa. In virtual colonoscopy, most mucosa voxels on the 

1 5 colon wall can be left intact in order to preserve more information. This can be very 
useful for three dimensional volume rendering. 

From the segmented colon wall volume, the inner surface, the outer 
surfece and the wall itself of the colon can be extracted and viewed as a virtual object. 
This provides a distinct advantage over conventional optical colonoscopy in that the 

20 exterior wall of the colon can be examined as well as the interior wall Furthermore, 
the colon wall and the colon lumen can be obtained separately from the segmentation. 

Because the colon is substantially evacuated prior to imaging, a 
conunonly encountered problem is that the colon lumen collapses in spots. While the 
mflation of the colon with compressed gas, such as air or COj, reduces the frequency 

25 of collapsed regions, such areas still occur. In performmg a virtual colonoscopy, it is 
desirable to automatically maintain a flight path through the collapsed regions and it is 
also desirable to use the scanned unage data to at least partially recreate the colon 
lumen in the coUapsed regions. Since the above described image segmentation 
methods effectively derive both the interior and exterior of the colon wall, this 

30 mformation can be used to enhance the generation of the fly path through the 
collapsed regions. 
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In extending the flight path through collapsed regions of the colon or 
expanding a collapsed region of the colon, the first step is to detect a collapsed region. 
Using the premise that the grayscale values of the image data &om around the outside 

r 

of the colon wall change much more dramatically than the greyscale values within the 
5 colon wall itself, as well as in other regions such as fat, muscle and other kinds of 
tissue, an entropy analysis can be used to detect areas of colon collapse. 

The degree of change in greyscale value, for example along the 
centerline, can be expressed and measured by an entropy value. To calculate an 
entropy value, voxels on the outer sur&ce of the colon wall are selected. Such points 
10 are identified from the above described image segmentation techniques, A 5x5x5 
cubic window can be applied to the pixels, centered on the pixel of interest Prior to 
calculating the entropy value, a smaller (3x3x3) window can be appUed to the pixels 
of interest in order to filter out noise firom the image data. The entropy value of a 
selected window about the pixel can then be determiaed by the equation: 

E=^C(i)ln{C{i)) 
i 

1 5 where E is the entropy and C(i) is the number of points in the window with the 

grayscale of i (i=0,l,2,. . 255). The calculated entropy values for each window are 
then compared against a predetermined threshold value. For regions of ak, the 
entropy values will be fairly low, when compared to regions of tissue. Therefore, 
along the centerline of the colon lumen, when the entropy values increase and exceed . 

20 the predetermiaed threshold value, a collapsed region is indicated. The exact value of 
the threshold is not critical and will depend iu part on the imaging protocol and 
particulars of the imaguig device. 

Once a collapsed region is detected, the previously determined 
centerline flight path can be extended tirough the region by piercing through the 

25 center of the collapse with a one voxel wide navigation line. 

In addition to automatically continuing the flight path of the virtual 
camera through the colon lumen, the region of colon collapse can be virtually opened 
usiag a physical modeling technique to recover some of the properties of the collapsed 
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region. In this technique, a model of the physical properties of the colon wall is 
developed. From this model, parameters of motion, mass density, damping density, 
stretching and bending coefficients are estimated for a Lagrange equation. Then, an 
expanding force model (i.e., gas or fluid, such as air, pumped into the colon) is 

5 formulated and applied in accordance with the elastic properties of the colon, as 

defined by the Lagrange equation, such that the collapsed region of the colon image is 
restored to its natural shape. 

To model the colon, a finite-element model can be applied to the 
collapsed or obstructed regions of the colon lumen. This can be performed by 

0 sampling the elements in a regular grid, such as an 8 voxel brick, and then applying 
traditional volume rendering techniques. Alternatively, an irregular volume 
representation approach, such as tetrahedrons can be applied to the collapsed regions. 

In applying the external force (air pumping) model to the colon model, 
the magnitude of the external force is first determined to properly separate the 

5 collapsed colon wall regions. A three dimensional growing model can be used to 
trace the internal and external colon wall surfaces in a parallel manner. The 
respective surfaces are marked from a starting point at the collapsed region to a 
growing source point, and the force model is applied to expand the surfaces in a like 
and natural manner. The region between the internal and external surfaces, i.e., the 

0 colon wall, are classified as sharing regions. The external repulsive force model is 
applied to these sharing regions to separate and expand the collapsed colon wall 
segments in a natural manner. 

To more clearly visualize the features of a virtual object, such as the 
colon, which is subjected to virtual examination, it is advantageous to provide a 

5 rendering of the various textures of the object. Such textures, which can be observed 
in the color images presented during optical colonoscopy, are often lost in the black 
and white, grey scale images provided by the CT image data. Thus a system and 
method for texture imaging during virtual examination is required. 

Figure 20 is a flow chart depicting a present method for generating 

0 virtual objects having a texture component The purpose of this method is to map 
textures obtained by optical colonoscopy images in the red-green-blue (RGB) color 
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space, as for example from the Visible Human, onto the gray scale monochrome CT 
image data used to generate virtual objects. The optical colonoscopy images axe 
acquired by conventional digital image acquisition techniques, such as by a digital 
"frame grabbef' 1429 which receives analog optical images from a camera, such as a 
5 video camera, and converts the image to digital data which can be provided to CPU 
1423 via interface port 143 1 (Figure 14). The first step in this process is to segment 
the CT image data (step 2010). The above described image segmentation techniques 
can be applied to choose intensity thresholds in the grey scale image to classify the CT 
image data into various tissue lypes, such as bone, colon wall tissue, air, and the like. 

10 In addition to performing image segmentation on the CT image data, 

the texture features of the optical image need to be extracted from the optical image 
data (step 2020). To do this, a gausian filter can be apphed to the optical image data 
followed by sub-sampling to decompose the data into a multiresolutional pyramid. A 
laplacian filter and steerable filter can also be applied to the multiresolutional pyramid 

IS to obtain oriented and non-oriented features of the data, While this method is 

elective at extracting and capturing the texture features, the implementation of this 
approach requires a large amount of memory and processing power. 

An alternative approach to extracting the texture features from the 
optical image is to utilize a wavelet transform. However, while wavelet 

20 transformations are generally computationally efficient, conventional wavelet 

transforms are limited in that they only capture features with orientations parallel to 
the axes and cannot be applied directiy to a region of interest. To overcome these 
limitations, a non-separable filter can be employed. For example, a lifting scheme can 
be employed to build filter banks for wavelets transform in any dimension using a two 

25 step, prediction and updating approach. Such filter banks can be synthesized by the 
Boor-Rom algorithm for multidimensional polynomial interpolation. 

After the textural features are extracted from the optical image data, 
models must be generated to describe these features (step 2030). This can be 
performed, for example, by using a non-parametric multi-scale statistical model which 

30 is based on estimating and manipulating the entropy of non-Gausian distributions 
attributable to the natural textures. 
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Once texture models are generated from the optical image data» texture 
matching must be performed to correlate these models to the segmented CT image 
data (step 2050). In regions of the CT image data where the texture is continuous, 
corresponding classes of texture are easily matched. However, in boundary regions 
5 between two or more texture regions, the process is more complex. Segmentation of 
the CT data around a boundary region often leads to data which is fuzzy, i.e., the 
results reflect a percentage of texture from each material or tissue and vary depending 
on the various weighting of each. The weighting percentage can be used to set the 
importance of matching criteria. 

10 In the case of the non-parametric multi-scale statistical model, the 

cross entropy or a Kullback-Leiber divergence algorithm can be used to measure the 
distribution of different textures in a boundary region. 

After texture matching, texture synthesis is performed on the CT image 
data (step 2050). This is done by ftising the textures from the optical image data in to 

15 the CT image data. For isotropic texture patterns, such as presented by bone, the 

texture can be sampled directly from the optical data to the segmented CT image data. 
For anisotropic texture regions, such as colon mucosa, a multiresolution sampling 
procedure is preferred. In this process, selective re-sampling for homogenous and 
heterogenous regions is employed. 

20 Alternatively, pseudocolor texture can be created directly from the CT 

data. For each voxel, multiple CT values, comprising a local area neighborhood, can 
be evaluated to determine a pseudocolor for the given voxel. For each voxel the local 
neighborhood consists of the voxels that are within some given distance of the center 
voxel. For example a 5x5x5 voxel cubic shaped region, or the double pyranodd which 

25 represents all voxels within 3 units measured by Manhattan distance. This vector of 
scalar values is then evaluated to map to a color to be displayed for this voxel during 
subsequent volume rendering. The evaluation of the local neighborhood vector of 
values can compute such things as local curvature, homo/heterogeneity, or other 
geometric or spatial functions. 



30 



Volume Rendering 
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In addition to image segmentation and texture mapping described 
above, volume rendering techniques can be lised in connection with virtual 
colonoscopy procedures to further enhance the fidelity of the resulting image. Figure 
21 illustrates a perspective volume ray-casting method which can be used for volume 

5 rendering in accordance with the present invention. From a selected virtual 
viewpoint, e.g., camera position, such as within the colon lumen, rays are cast 
through each of the proximate image pixels (step 2100). For each ray, the first 
sampling pomt is set as the current image pixel along the ray (step 2110). The 
distance (d) between the current sampling point and the nearest colon wall is then 

0 determined (step 2120). The current distance (d) is compared to a predetermined 
sampling interval (i) (step 2130). If the distance (d) is greater than the sampling 
interval (i) then no sampling occurs and the next sampling point along the ray is 
determined by jumping the distance d along the ray (step 21 40). If the distance is less 
than or equal to the sampling interval (i) then conventional sampling is performed on 

5 this point (step 2 1 50) and the next sampling point is selected in accordance with the 
sampling interval (i) (step 2160). For example, trilinear interpolation between the 
density values of 8 neighboring voxels can be performed to determine the new density 
value at the sampling point. 

The method of Figure 21 effectively accelerates ray-casting because a 

0 space leaping technique is used to quickly skip over empty space along the ray of the 
image plane to the colon wall. In this method, a distance firom a sample point to the 
nearest colon wall is determined along each ray. If the distance is larger than a 
predetemained sampling interval (i), a jump to the next sampling point along the ray is 
performed. Since the closest distance information is akeady available from the 

5 potential field which is used for virtual camera control, no additional distance coding 
calculations are required. In this case, neither surface rendering nor Z-huSex 
transform is required, which results in savuigs in preprocessing time and memory 
space* 

Altematively, a space leaping method can derive distance information 
0 for each ray from the Z-buffer of the corresponding surface rendering image. If the 
surface rendering image and volume rendering image will both be generated, this 
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* * 

approach provides minimal processing overhead burden as the Z-buffer infomiation is 
provided as a result of the surface rendering methods. Thus, this form of space leaping 
method only requires additional processing to perform a depth transformation from 
the image space domain to the world space domain. 
5 For those regions along the ray where the distance (d) was traversed in 

step 2140, the region along the ray corresponds to open space and can be assigned a 
value according to an open space transfer function. Typically, open space will have 
no contribution on the final pixel value. For each point where sampling takes place, 
one or more defined transfer fimctions can be assigned to map different ranges of 
10 sample values of the original volume data to different colors and opacities and 
possibly other displayable parameters. For example, four independent transfer . 
functions have been used to determine different material by mapping ranges of CT 
density values into specified colors of red, green, blue and opacity, each in the range 
of 0 to 255. 

15 Virtual Biopsy 

The above described techniques can also form the basis of a system for 
performing virtual electronic biopsy of a region being examined to effect a flexible 
and non-invasive biopsy. As noted above, volume rendering techniques use one or 
more defined transfer functions to map different ranges of sample values of the 

20 original volume data to different colors, opacities and other displayable parameters for 
navigation and viewing. During navigation, the selected transfer function generally 
assigns maximum opacity to the colon wall such that the outer surface is easily 
viewed. Once a suspicious area is detected during virtual examination, the physician 
can interactively change the transfer function assigned during the volume rendering 

25 procedure such that the outer surface bemg viewed becomes substantially transparent, 
allowing the region ioformation to be composited and thus the interior structure of the 
region to be viewed. Using a number of predetermined transfer functions, the 
suspicious area can be viewed at a number of different depths, with varying degrees of 
opacity assigned throughout the process. 
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Polyp Detectioii 

The present system and methods can be used to perform automated 
polyp detection. With reference to Figure 13, polyps 1305, which occur, for example, 
within the colon, generally take the form of small convex hill-like structures 
5 extending from the colon wall 1301. This geometry is distinct from the fold of the 
colon wall. Thus, a differential geometry model can be used to detect such polyps on 
the colon wall. 

The siuface of the colon lumen can be represented as a continuously 
second diff^entiable surface in three dimensional Euclidean space, such as by using a 

10 C-2 smoothness surfece model. Such a model is described in "Modem Geometry 
Methods and Applications" by B.A. Dubrovin et al, published by Springer-Verlag 
1994, which is hereby mcoiporated by reference in its entirety. In this model, each 
voxel on the surface of the colon has an associated geometrical feature which has a 
Gauss curvature, referred to as Gauss curvature fields. A convex hill on the surface, 

1 5 which may be indicative of a polyp, possesses a unique local feature in the Gauss 
curvature fields. Accordingly, by searching the Gauss curvature fields for specific 
local features, polyps can be detected. Once detected, the suspected polyps can be 

■ 

highUghted and thus brought to the attention of the physician where the physician can 
measure the suspected polyp and use the above described virtual biopsy methods to 
20 further investigate the suspicious region. 



Central Fly-Path Generation 

In the case of virtual colonoscopy, determining a proper navigation 
line, or fly-path, through the colon lumen is an important aspect of the described 
systems and methods. While certain techniques for determining the fly-path of the 

25 virtual camera model were discussed with respect to Figures 4-8, Figure 22 illustrates 
an alternate method of generating the central fly-path through the colon lumen. After 
the colon wall is identified, such as by the image segmentation methods described 
herein, a volume shrinking algorithm can be employed to emphasize the trend of the 
colon lumen and reduce subsequent searching time within the lumen volume (step 

30 2310). 
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Figure 23 further illustrates the steps of an exemplary volume 
shrinking algorithm, which is based on a multiresolution analysis model. In this 
procedure, the three dimensional volume is represented by a stack of binary images 
which have the same matrix size (step 23 10). Collectively, these images form a 
5 binaiy data set. A discrete wavelet transformation can be applied to the binary data 
set which results in a number of sub-data sets representing different time-frequency 
components of the binary data set (step 2320). For example, the discrete wavelet 
transformation may yield eight (8) sub-data sets. The sub-data sets are compared 
against predetermined threshold values such that the lowest frequency component is 

10 identified (step 2330). This component forms the binaiy data set for subsequent 

discrete wavelet transfom[iation and thresholding steps, which are recursively applied 
in a multi-resolution structure (step 2340). In the case of virtual colonoscopy, the 
discrete wavelet transformation and associated thresholding can be applied three times 
recursively on the subsequent sub-dataset that represents the lowest frequency 

1 5 component (a 3 -level multi-resolution decomposition). 

Returning to Figure 22, from the reduced colon volume model, a 
distance map technique can be employed to generate a Tninimnm distance path 
between the two ends of the colon, e.g., from the rectum to the cecum (step 2215). 
The resulting path preserves the global trend information of the colon lumen, but 

20 ignores the trends exhibited by local folds. Control points within the global colon can 
then be determined by mapping the minimum distance path back to the original data 
space (Step 2220). For example, in the case of a 3-level multi-resolution 
decomposition, the reduced volume is three times smaller than the original volume 
and an affine transformation, which is well known, can be used to map the reduced 

25 volume model exactly back to the original scale volume. The minimum distance path 
of the reduced value can also be mapped back into the original scale volume as a 
series of points, which can be used as the control points within the colon. 

The preferred fly path is one which is on the centerline of the colon 
lumeiL However, the initial control points may not be exactiy located in the center of 

30 the colon Imnen. Thus, the initial control points can be centered, such as by the use of 
a bi-section plane algorithm (step 2230). For example, at each selected control point. 
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a bi-section plane can be defined as a plane nonnal to the trend direction and cutting 
across the colon lumen. A centralization algorithm, such as a maximum disk 
algorithm, can then be performed on each bi-section plane. Such an algorithm is 
discussed in the article "On the Generation of Skeletons from Discrete Euclidean 
5 . Distance Maps" by Ge et al., IEEE Transactions on PAM, Vol. 18, pp. 1055-1066, 
1996 which is hereby incorporated by reference. 

Once the control points are centralized, the flight path can be 
determined by interpolating a line connecting these points (step 2240). In the case of 
virtual colonoscopy, it is desirable that the interpolated flight path take the form of a 

1 0 smooth curve which is substantially centered within the colon lumen. A constrained 
cubic B-spline interpolation algorithm based on Serret-Frenet Theorem in differential 
geometry theory can be used to establish a suitable smooth curved flight path, such as 
is described in '"Numerical Recipes in C: The Art of Scientific Computing," by Press 
et al.. Second Edition, Cambridge University Press, 1992. 

1 5 The pictorial representation of a segmented colon lumen in Figures 24 

and the flow chart of Figure 25 set forth yet another alternate fly-path generation 
method in accordance with the present invention. In this alternate method, the 
representation of the colon lumen 2400 is first partitioned into a number of segments 
2402 a-g along the length of the lumen 2400 (step 2500). From within each segment 

20 2402 a representative point is selected 2404 a-g (step 2520). Each representative 

point 2404 a-g is then centered with respect to the colon wall (step 2530), such as by 
the use of a physically-based deformable model which is used to push the points to the 
center of the respective segment. After the representative points are centered, the 
points are sequentially joined to establish the center-line fly-path for the virtual 

25 camera model (step 2540). If tiie segments are sufficiently small in length, the 
centered points can be connected with straight Ime segments 2406 a-f. However, 
when linear curve fitting techniques are applied to join the centered pomts, a 
smoother, continuous flight path is established. 

Each of the foregoing methods can be implemented using a system as 

30 illustrated in Figure 14, with appropriate software being provided to control the 
operation of CPU 1409 and CPU 1423. 



wo 02/029764 PCT/USO 1/30704 

45 

An alternate hardware embodiment, suitable for deployment on a 
personal computer, is illustrated in Figure 26. The system includes a processor 2600 
which should take the form of a high speed, multitasking processor such as a Pentium 
in processor operating at a clock speed in excess of 400 MHZ. The processor 2600 is 
coupled to a conventional bus structure 2620 which provides for high speed parallel 
data transfer. Also coupled to the bus stmcture 2620 are main memory 2630, a 
graphics board 2640, and a volume rendering board 2650. The graphics board 2640 is 
preferably one which can perform texture mapping, such as the Diamond Viper v770 
Ultra manufactured by Diamond Multimedia Systems. The volume rendering board 
2650 can take the form of the VolumePro board from Mitsubishi Electric, which is 
based on U.S. Patent Nos. 5,760,781 and 5,847,71 1, which are hereby incorporated by 
reference. A display device 2645, such as a conventional SVGA or RGB monitor, is 
operatively coupled to the graphics board 2640 for displaying the image data. A 
scanner interface board 2660 is also provided for receiving data from an imaging 
scanner, such as an MRI or CT scanner, and transmitting such data to the bus structure 
2620. The scanner interface board 2660 may be an application specific interface 
product for a selected imaging scanner or can take the form of a general purpose 
input/output card. The PC based system 2600 will generally include an I/O interface 
2670 for coupling I/O devices 2680, such as a keyboard, digital pointer (e.g., mouse) 
and the like to the processor 2620. Alternatively, the I/O interface can be coupled to 
the processor 2620 via the bus 2620. 

In the case of three dimensional imaging, including texture synthesis 
and volume rendering, nimierous data handling and processing operations are 
required. For large datasets, such as those represented by the colon lumen and its 
surrounding area, such processing can be very time consuming and memory intense. 
However, using the topology of Figure 26 in accordance with the processing method 
illustrated in the flow chart of Figure 27, such operations can be performed on a 
relatively low cost personal computer (PC). Imaging data is received by the processor 
2620 and stored in main memory 2630 via the scanner interface board 2660 and bus 
structure 2620. This image data (pixels) is converted into a volume element (voxel) 
representation (step 2710). The volume representation, which is stored in main 
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memory 2630, is partitioned, for example into slices, such as along a major volume 
axis or other portions of the region being imaged (step 2720). The volume partitions 
are then transferred to the volume rendering board and temporarily stored in volume 
rendering memory 2655 for volume rendering operations (step 2730). The use of 

5 locally resident volume rendering memory 2655 provides for enhanced speed in 
volume rendering as data need not be exchanged over the bus 2620 during rendering 
of each slice of the total volume. Once volume rendering is complete for the slice, the 
rendered data is transferred back to main memory 2630 or the graphics board 2640 in 
a sequential buffer (step 2740). After all slices of interest have been subjected to 

0 rendering, the contents of the sequential buffer are processed by the graphics board 
2640 for display on the display unit 2645 (step 2750). 



Multi-scan Based Virtual Examination 

The techniques discussed above generally perform virtual imaging 
based on a dataset acquired &om a single magnetic resonance imaging (MRI) or 
computed tomography (CT) scan. However, the techniques discussed above are also 
useful for performing virtual examination of a region using multiple scans of a region. 
By using multiple scans of a region, improved imaging of regions of pathology can be 
achieved and motion artifacts can be reduced. One such appUcation of interest is in 
performing virtual cystoscopy to screen a patient for possible polyps or cancer of the 
bladder. 

Figure 28 is a flow chart which illustrates a method of employing 
multiple MRI scans to perform virtual examination of an object, such as virtual 
cystoscopy. Unlike CT images, where the bladder wall can be difficult to distinguish 
from urine, in MRI images, urine can be used as a natural contrast agent to delineate 
the inner bladder wall. To this end, a pre-image scan protocol is employed (step 2805). 
Approximately V2 hour prior to the first of four MRI scans, the patient is requested to 
empty the bladder and then consume one cup of water. After approximately 14 hour, 
the patient is subjected to the first of four MRI scans of the bladder region (step 2810). 
The first scan, with the bladder full and distended, follows protocol for Tl -weighted 
transverse imaging. For example, when using the Picker scanner referenced above, a 
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KJELL FASTER protocol using a 256x256 matrix size, a 38 cm field of view (FOV), a 
1 .5 mm slice thickness (no gap), a 3 ms TE, a 9 ms TR, a 30 degree flip angle and one 
scan average can be used. Of course, these parameters tend to be scanner specific and 
various changes in the parameters can be used with acceptable results. 

5 With the bladder still fiiU, the patient is subjected to a second MRI scan, 

scan 2 (step 2815). The second scan follows a protocol for Tl-weighted coronal 
imaging, such as the Picker KJELL FASTER protocol with a 256x256 matrix size, a 
38 cm field of view (FOV), a 1 .5 mm slice thickness (no gap), a 3 ms TE, a 9 ms TR, a 
30 degree flip angle and a two-scan average. 

0 The two image scans described above are taken along orthongal axes 

with respect to one another. The advantage of this is that regions of significant motion 
artifacts in one scan, generally correspond to regions of minimal motion artifacts in the 
orthogonal scan. Accordingly, by taking a first scan in the transverse direction and a 
second scan in the coronal direction, the image scans can be registered and motion 

5 artifacts in the data set can be identified and compensated for. 

After the scan 2, the patient is asked to relieve the bladder and is then 
subjected to two additional MRI scans. The third scan (step 2820) follows the same 
imaging protocol as the first scan (transverse imaging). The fourth scan (step 2825) 
follows the same imaging protocol as the second scan (coronal imaging). 

0 The image scans can be acquired using a Picker 1 .5 T Edge whole-body 

scanner. Although a T2 imaging protocol can be used, a Tl imaging protocol is 
preferred for virtual cystoscopy because this protocol provides improved delineation 
between fat and urine and offers a shorter acquisition period. Alternatively, the image 
scans can take the form of computed tomography or ultrasound imaging scans using 

5 suitable contrast agents and protocols for these imaging techniques. 

During the first two scans (scan 1 and scan 2), the bladder is distended 
and the bladder wall is relatively thin. In this case, physiologically altered locations, 
such as tumors, may thin at a different rate as compared to the unaltered bladder wall 
and may become more apparent under these conditions. During the third and fourth 

0 scans, the bladder is substantially empty and the bladder wall is thicker. With a thicker 
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wall, a more pronounced image contrast may result between normal tissue of the 
bladder wall and that of physiologically altered tissue. 

After the four scans are acquired, the four corresponding datasets can 
then be individually processed. Initially, each scan data set is preferably subjected to 

S image segmentation, as discussed above (step 2830), such as in connection with Figure 
15. During image segmentation, the voxels of the four datasets are classified into a 
number of categories, such as bladder wall, urine, fat, boundary, etc. The classification 
is based on the local intensity vectors of the voxels. Once the voxels are classified, the 
interior of the bladder lumen can be identified using a region growing algorithm 

0 beginning with a seed voxel selected firom within the bladder volume, such as by 
selecting an air voxel or urine voxel. 

Prior to clinical analysis of the segmented volume data sets, registration 
of the four data sets to a common coordinate system is performed (step 2835). Because 
the shape of the bladder varies firom scan to scan, an exact voxel-voxel registration is 

5 not of practical value. Instead a flexible registration process is preferred. In the 

present flexible registration process, for each volume of interest (volume rendered for 
each corresponding scan) the center of the volume is determined, such as by averaging 
the three coordinates of all the voxels in the volimie. 

A Cartesian coordinate system can then be constmcted with the origin 

0 of the system located at the center point of the volume. The axes of the system can 
then be oriented in a number of ways. A suitable selection of orientation corresponds 
to the orientation of the natural human body, e.g., with the Z-Axis running along the 
height of the body (e.g., fi:om toe to head) the Y-axis oriented firom back to firont and 
the X-axis running laterally ( e.g., firom left to right). The units of length in this 

5 coordinate system can be convenientiy set to an arbitrary unit of one voxel length, the 
absolute magnitude of which will vary based on acquisition properties for the MRI 
scans. So long as the same pixel spacing is used in all scans to acquire all four data 
sets, this will result in a uniform value for each of the four data sets. 

After registration, the images firom the four data sets can be viewed 

0 individually or simultaneously (step 2845). An exemplary display window is 

illustrated in Figures 29 and 30. Referring to Figure 29, the display is partitioned into 
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four sub windows 2905, 2910, 2915, 2920 which correspond to scan 1, scan 2, scan 3 
and scan 4, respectively. A control panel section 2925 can also be provided on a 
portion of the display to establish a graphical user interface (GUT) to offer display and 
navigation functions to the user. As an operator navigates in one of the unage sub 
5 windows, such as magnifying the view, the correspondii^ operation preferably takes 
place in the other sub window views as well. A user can also select one of the views 
for expansion to a single window display. 

To reduce the amount of data which is simultaneously processed, the 
data sets can be partitioned, such as into 8 parts or octants (step 2840). This can be 

1 0 performed in a number of ways. For example, with reference to the Cartesian 

coordinate system illustrated in Figure 29, the data can be portioned into the eight 
regions of the coordinate system: (1) X, Y, Z; (2) X, -Y, Z; (3) X, Y, -Z; (4) X, -Y, -Z; 
(5) -X, Y, Z; (6) -X, -Y, Z; (7) -X, Y, -Z; and (8) -X, -Y, -Z. 

Figure 29 illustrates four views of the outside of the bladder lumen 

1 5 taken from each of the four scans. Figure 3 0 illustrates fours views of a portion of the 
interior of the bladder lumen also taken from each of the four scans. 

Multi-Resoltttion Ima^ng and Virtual Laryngoscopy 

The systems and methods described herein can be adapted and applied 

20 to perform multiresolution imaging which is well suited for virtual laryngoscopy. 
Figure 3 1 is a flow chart illustrating a method for performing virtual laryngoscopy. 
First, an imaging scan of the region of a patient's larynx is acquired (step 3105). This 
can be performed using computed tomography (CT) or magnetic resonance imaging 
(MRI) techniques. However, because the CT scan hi this region offers significantly 

25 faster acquisition time (30 seconds versus over 7 minutes for MRI) and higher 

resolution (0.3mm cubic voxel compared to 1mm cubic voxel for MRI), the CT scan is 
preferred. To acquire the CT scan data a GE/CTI spiral scan CT scanner can be used. 
A suitable scan protocol is 120 keV, 200 ma, 512x512 matrix size, 15 cm field-of-view 
and 3mm/2.0:l pitch. The scan is completed in approximately 30 seconds and results 

30 m 351 image slices of 0.3mm thickness and results in 0.3mm cubic voxels. 
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Image segmentation can be used to classify voxels into a number of 
categories (step 3110). In this operation, a modified self-adaptive on-line vector 
quantization (SOVQ) algorithm can be used. In such a case, the algorithm analyzes 
each voxel vnth respect to neighbors of up to the third order to determine local density 
5 features. Each voxel in the acquired dataset has an associated local density vector. By 
transforming the local density vectors using the Karhunen-Loeve (K-L) transform, 
feature vectors for the voxels in the volume image can be obtained. Based upon the 
feature vectors, the voxels can be classified and labeled. Voxel classification is 
dependent in part on the choice of a local voxel density vector and one preset 

1 0 parameter, referred to as the maximum cluster number (MCN). The MCN sets the 
number of voxel classifications that will be applied to the dataset. In the case of the 
CT images, the human eye can discern four (4) distinguishable tissue/material types. 
An MCN value of 5 is suitable in this case. For an MRI image, the human eye can 
differentiate among 6 different tissue types, and an MCN value of 7 can be used. 

IS As part of the image segmentation process, an expanded data set is 

generated by interpolation between the measured data points. For example, prior to 
employing the SOVQ algorithm, a first order Lagrange interpolation can be applied to 
each slice in the dataset. This expands the 256x256 matrix size of the original slices of 
the data settoa512x512 matrix size. In addition, inter-slice interpolation can be ' 

20 perforaied to further expand the dataset between actual slices. The interpolated dataset 
is referred to as the enlarged dataset. In addition to generating an enlarged dataset, the 
interpolation process also suppresses noise and reduces the partial-volume effect, as 
the interpolation process has a low-pass filtering effect on the data. 

Using a two dimensional viewing tool, a seed voxel can be selected 

25 within the larynx lumen and a growing algorithm applied to extract the larynx volume 
fi-om the dataset (step 3115). In those regions of the larynx where there may be several 
uncormected volume regions, multiple seed points can be selected. 

With the larynx volume identified and the voxels of the regions 
classified through image segmentation, the next task is to manage the data in a manner 

30 which allows efficient navigation and viewing of the virtual larynx. In this case, a 

level-of-detail (LOD) approach is adopted and modified for use in the present method. 



wo 02/029764 PCT/USO 1/30 704 

51 

In this LOD method, a reduced dataset is generated from the enlarged data set For 
example, the 512x512x256 enlarged dataset can be reduced to a 64 x 64 x 32 reduced 
volume dataset using a multi-resolution decomposition with three levels of 
tfiresholding (step 3 120). Next, polygons used to render the volume unages in both the 
5 enlarged and reduced volume datasets can be extracted. A traditional Marching Cubes 
method can be used to extract polygons to fit the surface of the larynx. 

One problem encountered in the prior art is managing the large number 
of polygons required to generate the three dimensional image for the enlarged dataset 
This problem is solved in the present method by organizing the enlarged dataset in a 

10 Binary Space Partitioning (BSP) tree data structure (step 3130). The original image 
volume is selected as the root of the tree. The space is then partitioned into two 
subspaces contaiiiing an approximately equal number of polygons. This subdivision 
process is iteratively repeated until the number of polygons in each resulting subspace 
is below a threshold value. The threshold value can vary based on system performance 

1 5 and application requirements. The last resulting subspaces are referred to as leaf nodes 
of the tree. Once the subdivision process is complete, all of the voxels of the expanded 
dataset are stored in the leaf nodes of the BSP tree. 

During navigation or viewing, polygon culling can be appUed by first 
removing those leaf nodes that are completely outside the field-of-view from current 

20 processing operations. The remaining polygons are recalled from the BSP tree, 
ordered and rendered in those spaces which were not culled. Thus, the BSP tree 
provides an effective tool for selecting a relevant portion of the dataset for a particular 
navigation or display operation. 

The enlarged and reduced datasets are cooperatively used in a two level 

25 LOD rendering mode. If a user is interacting with the object (step 3 1 35), such as 

rotating, shifting or effecting other changes in the field of view, the polygons from the 
reduced dataset (64-sized) are rendered (step 3140). Because of the significantly lower 
number of polygons involved, interaction with the reduced dataset volume can be 
performed faster and with less processing overhead. The tradeoff for the increased 

30 speed is reduced image resolution. If there is no interaction from the user after a 
predetermined time period, the polygons of the enlarged dataset (5 1 2-sized) are 
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selected from the BSP tree and are rendered to provide a high resolution image of the 
current field of view (step 3145). 

Virtual Angiography 

The techniques for virtual imaging and navigation can also be adapted 
5 and appUed to virtual angiography. This technique can be used for detection and 

measurement of various abnormalities and disease of the chrculatory system. 

One such application of virtual angiography is the detection of 

abdominal aortic aneurysms, which generally start as small enlargements of the aortic 
• vessel and exhibit a greater risk to rupture with increasing size of the aneurysm. 
0 Previously, the only effective method of treatment was open surgery, placing a graft 

within the aorta at the level of the aneurysm. However, this procedure has a high 

degree of associated morbidity and mortality. Recently developed per cutaneous 

placed aortic stent graft techniques have a significantly lower complication rate. 

Virtual angiography is an effective method to help plan these less invasive procedures 
5 and can also be an efifective tool for detecting the presence of an aneurysm and tracking 

the growth of an aneurysm to determine if and when surgery is indicated. 

Figure 32 in a flow chart which provides an overview of the present 

virtual angiography method. In performing a virtual angiography, an image scan of the 

vessel, such as the aorta must be acquired (step 3205). Various imaging techniques 
0 can be used, such as Computed Tomography (CT), Magnetic Resonance Imaging 

(MRT) and ultrasoxmd. However, an aortic CT scan is generally preferred because of 

the contrast between blood, soft tissue and calciimi deposits which results in the CT 

image. 

Once an unage scan data set is acquired, image segmentation techniques 
5 are then applied to the data set to classify the voxels of the dataset into a number of 
categories (step 3210). The image segmentation techniques described above, such as 
in connection with Figure 15, are generally applicable. In this case, the various feature 
vector values of the voxels will be grouped according to categories such as blood, soft 
tissue and calcium deposits. Using a blood voxel as a seed, a region growing algorithm 
0 can be used to determine the volume and extent of the aortic lumen. 
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In the CT image, an aneuiysm has image features which closely 
resemble the neighboring soft tissue. As a result, the fidl contour of the aneurysm can 
be difficult to establish. However, regions with calcium deposits offer significant 
contrast on the CT scan and can be used to identify portions of the aneurysm, such as 
5 the endpoints of the aneurysm on the vessel wall (step 3215). 

After a portion of an aneurysm is detected, one or more closing siirfaces 
can be generated to define an estimation of the aneurysm's contour (step 3220), A 
convex closing surface can be established using a non-uniform, non-rational B-spline 
to generate a surface which runs through or near the points of the aneurysm which were 
10 identified. 

After the closing surface is generated, the volume of the aneurysm can 
be estimated (step 3225). One method for estimating the volume is to count the 
number of voxels which are enclosed by the estimated closing surface. In addition, 
within the volume of the aneurysm, the centerline along the direction of blood flow can 

15 be determined by using a distance transform technique. Continuous local coordinate 
systems can then be established along this centerline and the diameter of the aneurysm 
determined. Virtual navigation can take place along this centerline, in a manner 
consistent with that described above for navigating through a lumen, such as the colon. 

Referring to Figures 33A-C, the described method of virtual 

20 ^ angiography can be used to assist in the generation and placement of a stent graft to 
bypass an abdominal aortic aneurysm. Figure 33 A illustrates a simplified diagram of 
an abdominal aortic aneurysm located below the renal arteries and above the 
bifiircation of the aorta. Because of variations from patient to patient in the specific 
anatomy of the aorta and the size and location of an abdominal aortic aneurysm therein, 

25 when a stent graft is to be used to bypass an aneurysm, the graft must be designed and 
built to specifically fit the particular aortic segment. As illustrated in Figure 33B, this 
can require identifying the length of the required graft, the diameter at the points of 
interface on each end of the bypassed region, the angles of interface, among other 
variables. If the aneurysm is located near an arterial branch, the size and angles of the 

30 bifijrcated ends of a bifiircated stent graft must also be determined, as illustrated in 
Figures 33B and 33C. 
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To date> such measurements have been performed through invasive 
calibrated angiograms using a catheter inserted into the aorta from an insertion made at 
the level of the groin region, rapid injection of a large amount of iodinated contrast and 
rapid radiographic imaging. This technique can be supplemented and perhaps 
5 supplanted using the present virtual angiography techniques, which can resolve such 
distances and angles using virtual navigation using centerlines constructed through the 
branches of the aortic lumen. In addition, the virtual angioscopy can be used to 
perforai a virtual biopsy of the region where a stent graft may be inserted. This allows 
the operator to view beneath the arterial surface and examine the region for thrombus 
10 deposits, calcification or other factors which would contra-indicate the use of a stent 
graft procedure. 

Another application of virtual angiography is the imaging, examination 
and navigation through the carotid arteries which supply blood flow to the brain. The 
techniques described herein with respect to virtual endoscopy are fully applicable in 

15 the case of blood vessels. For example, the vessels of interest are extracted fiom the 
acquired image data using image-segmentation techniques. Next, a navigation flight 
path can be established through the vessel(s). Preferably, potential fields are built up 
within the vessel for use in navigation. As with other organs, such as the colon, a 
volume-rendered model of the vessels of interest can be generated. Using the flight 

20 path and potential fields to navigate through the interior of the volume rendered blood 
vessel lumen, abnormalities such as vessel narrowing and plaque build up can be 
observed. In addition, the techniques discussed regarding virtual biopsy can be applied 
in this context to evaluate vessel wall and characterize build up on the wall surface, 
such as plaque. 



25 

Tree Branch Searching for Virtual Endoscopy 

Path planning for virtual navigation through a hollow organ or object is 
an important task. Various techniques have been discussed, such as fly-path 
generation, to achieve this goal. As the geometry of the object being studied becomes 
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more complex, such as presenting a multi-branch structure^ the task of path planning 
becomes even more complex. It is desirable to determine not only the center line of a 
primary lumen, but also to identify and locate any branches extending from the primary 
lumen. Common examples of organs having a complex branch structure include the 
5 main airway and lungs, the cardiovascular system and, because of the presence of 
haustral folds, the colon. Each organ or object generally presents specific challenges 
for defining a path, or skeleton, for the object. However, a generalized technique for 
generating such a skeleton is illustrated in the flow chart of Figure 34. 

Referring to Figure 34, an imaging scan of the region of interest, such as 

10 a computed tomography (CT) or Magnetic Resonance Imaging (MRI) scan, is acquired 
(step 3405). As discussed above, the imaging scan is transformed into a three 
dimensional volume of the region by stacking the binary images of the imaging scan 
and defining three dimensional volume units, or voxels, from these stacked images 
(step 3410). Depending on the volume and complexity of the region of interest, it may 

15 be desirable to reduce to size of the dataset of the three dimensional volume prior to 
generating the skeleton. To this end, a multiresolution data reduction process, which is 
discussed in more detail below, can be used (step 3415). 

The skeleton is a subset of the three dimensional volume. Preferably, 
the skeleton has the following attributes: 1) It preserves the homotopy of the tree; 2) it 

20 is 26-<:onnected; 3) it is one voxel thick; 4) it approximates the central axes of the 
branches; and 5) it is relatively smooth. The degree of homotopy is somewhat 
application specific. For example, m generating a skeleton of the colon lumen, the 
skeleton will generally be a single path from end to end, despite the presence of 
numerous haustral folds which can be several voxels deep. However, in the 

25 cardiovascular system and pulmonary system, a small offshoot from the root which is 
several voxels deep can represent a legitimate branch in the system. 

Returning to Figure 34, in the volume of interest, a root voxel is 
identified in the volume (step 3420). In performing virtual endoscopy, this can be 
performed manually based on an understanding of the geometry of the structure being 

30 evaluated. 



wo 02/029764 



PCT/USOl/30704 



56 

A distance map can then be generated to identify the branches in the 
tree and the distances between the endpoints of the branches and the root voxel (step 
3425). A presunaption applied in this method is that there exists one unique endpoint 
on each branch which exhibits the longest distance to the root of the tree. Figure 35 is 
5 a schematic diagram illustrating a 3x3x3 cubic voxel arrangement which is referred to 
as a 26-connected voxel cubic distance plate. In the center of this arrangement is a 
seed voxel 3505, which is assigned a distance weight of zero. Around the seed voxel 
3505 are 26 connected neighbor voxels which are assigned distance weights based on 
the respective Euclidean distance between the respective neighbor voxel and the seed. 
10 In a cubic arrangement the Euclidian distance can assume a normalized value of 

which is approximately equal to 1, 1.4 and 1.7. To simplify processing, 

the voxels can be assigned integer value weights of 10, 14, and 17 to approximate the 

relative Euclidian distances. 

Figure 36 is a pseudo-code representation of an algorithm for 
1 5 determining the distance map from a voxel in the volume to the root using the 

Euclidian weighted distances of the 26-comiected cubic distance plate of Figure 35. 

From the generated distance map, branches are identified and the endpoints of the 

branches are determined (step 3430). 

Referring to Figure 36, the root of the volume is labeled with the integer 
20 value 0. A processing queue is then formed with the voxels in the volume. The 

voxels are then labeled in a first-in, first out maimer by adding the Euclidian distances 

between the voxel at the top of the queue and the root voxel. This process is repeated 

until all of the voxels in the volume are assigned a value in the distance map. 

Because the labelmg of voxels in the distance map will depend, in part, 
25 on the queuing order, the resulting distance map does not provide a unique solution. 

However, regardless of the queuing order, there is always at least one farthest point for 

each branch. In addition, for each voxel, other than the root voxel, there is always at 

least one 26-connected neighbor in the volume which has a shorter distance to the root. 

Thus, the endpoints of the branches are readily detectable by searching the distance 
30 map for local maximum distances (local maxima) (step 3430). The term local maxima 

is a relative terai. In evaluating the volume for local maxima, the volume should be 
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partitLoned into various subspaces which are appropriate to the object being evaluated. 
The expected feature size, branch length, branch diameter, etc. are generally considered 
in detennining the subspace partitions. 

Once the endpoints of the branches are determined, the shortest path 
5 from the endpoint to the root voxel is detennined (step 3435). The shortest paths from 
the endpoints to the root define the basic structure of the branches of the tree and 
approxinotate the centeriine of the branches. This is referred to as the rough skeleton of 
the volume. The shortest paths are preferably generated from the braQches at farthest 
end of the tree and begin from the end of those branches. From the most remote 

10 branch endpoint, the first voxel is selected and its 26-connected neighbors are analyzed 
to detennine which voxel is in the minimal distance path from endpoint to root. This 
process is repeated until a selected voxel meets the root. This results in a one-voxel 
wide path from the farthest end to the root Searching for the shortest path for other 
branches is similar. However, for subsequent branches, the selection process can 

1 5 terminate when the current path reaches a previously assigned path (e.g., the path need 
not lead all the way to the root). The collection of all of the interconnected shortest 
paths is the resulting rough skeleton of the object. 

Depending on the application of the resulting rough skeleton, it may be 
desirable to refine the rough skeleton (step 3440). One step of refining the skeleton is 

20 to centralize the skeleton within the branches. Centralization preferably takes place 
branch by branch from the longest branch to the shortest Starting with the longest 
branch, a uniform interval is selected, generally in the range of 4-8 voxels, along the 
branch. For each interval, the tangent direction of the voxel on the rough skeleton is 
calculated and a plane crossing the voxel perpendicular to the tangent direction is 

25 determined. A two dimensional area defined by the intersection of the plane and the 
volume is created and the center of this intersection can be computed using the known 
maximum disk technique. The centers of intersection can then be connected using a 
bi-cubic, B-spline interpolation or other curve fitting method. For the remaining 
branches, the endpoint which meets another branch or the root must first be adjusted to 

30 match the position of the previously centered skeleton branch. Then, centralization can 
proceed in the same manner as described for the longest branch. 
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Referring back to step 3415, when a large dataset is involved, it may be 
required, or at least desirable, to reduce the size of the dataset to speed up processing 
and reduce processing cost. Noting that the tree structure can be preserved within a 
range of scales, the large volume can be shrunk to a smaller scale space for structure 
5 analysis. 

A shrinking method based on multiresolution analysis theory can be 
used. The input data is the stack of binary images of the same size which can be 
obtained from the segmentation results of the CT or MRI scan. The x-direction is 
taken along the slice image width, the y-direction is along the slice image height, and 

» 

0 the z-direction is along the direction of slice by slice. The foreground voxels in the 
tree volume are set to value of 128 (maximum) and the background voxels are set to 
value 0 (minimum). A Daubechies' bi-orthogonal wavelet transform with all rational 
coefiBcients is employed. This one-dimensional (ID) discrete wavelet transformation 
(D WT) is first applied along to the x-direction row by row. From application of the 

5 DWT only the lower frequency components are retained and packed. The computation 
is preferably implemented in floating points. Noting that the DWT is applied to the 
binary signal, there are two kinds of nonzero coefficients which result in the lower 
frequency component. The first is of value 128 and this kind of coeffrcients are located 
in the interior of the volume. The second is of a value not equal to 128 and these 

;0 coefficients locate the boundary of the volume. 

The coefficients of the second kind are compared against a 
predetermined threshold value. If its absolute value is larger than a pre-set threshold 
Tl, the value of the coefficient is set to 128; otherwise, it is set to 0.. This results in a 
stack of binary images with a row size of half of the original dataset. The same DWT 

5 is then applied to the resulting dataset along the y-direction colunm by column, where 
the similar thresholding is employed to the lower frequency components. The result is 
again a stack of binary images, but now with both half row and colxmm size as 
compared to the original dataset. Finally, the DWT is applied to the last result along 
the z-direction and the lower frequency components are retained. This step completes 

1 0 the first level decomposition. 
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The resulting dataset of the first level decomposition is of half size in all 
three directions as compared to the original dataset. If the shrinking procedure stops at 
this level, the finial thresholding is applied. It revalues those coefficients of nonzero 
and non-128 value. If the absolute value of this kind of coefficient is larger than a pre- 
5 set threshold T2, it will be revalued as 128; otherwise, it is revalued as 0. If further 
shrinking is needed, the same thresholdmg algorithm is applied with the threshold Tl. 
Further shrinking proceeds as previously described, but is applied to the dataset shrunk 
at the last previous level. The decomposition procedure can be recxirsively applied 
until the resulting volume meets the desired reduced data volume. In virtual 

10 endoscopy, the slice unages are of 512X512 pixel size. The maximum decomposition 
level is usually three, resulting in a 64x64 reduced pixel size. 

The volume is isotropically shrank in all directions with the presented 
method. The two pre-set thresholds, Tl and T2, are used to control the degree of 
shrinking. If the volume is significantly over shrunk, connectivity may be lost in the 

1 5 reduced volume. If it is over shrunk to a leaser degree, two separate branches may 
merge into one branch in the reduced volume dataset The larger the two threshold 
values, the thinner the reduced volume is. The range of those two thresholds is [0, r x 
128], where 0<Kl. Preferably, the range for virtual endoscopy is re (0.08, 0.28) for 
Tl and re (0.7, 0.98) for T2. The exact determination is dependant on the feature size 

20 of the particular application and is selected to achieve reduction while retaining the 
fidelity of the structure information in the shnmk volume. 

After shrinking the original volume, the tree branch searching procedure 
can be applied to the smaller volume (steps 3420-3440). The resultant skeleton can be 
mapped back into the original scale space. When scaled to the origmal space, the 

25 image of the smaller scale skeleton no longer remain a connected path in the original 
scale space. These voxels in the image act as control pomts for the final skeleton. The 
control points are centralized using the algorithm as described previously, and then, 
they are interpolated to form the final skeleton of the object. 

Computer Assisted Diagnosis 
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The virtual examination techniques described herein lend themselves to 
applications for the computer assisted diagnosis (CAD) of various conditions. For 
example, as described above, by examining the geometry of an organs tissue for local 
Gausian curvatures, regions with abnormal geometry, such as polyps inside a colon 
S lumen, can automatically be identified. This technique can be generalized and used in 
conjunction with texture features to provide CAD functionality for a number of 
applications. 

For example, using the multi-scan imaging of the bladder described 
above, automated detection or tumors in the bladder wall can be performed. In this 

0 case, the degree of tumor invasion within the bladder wall is generally used to define 
the stage of bladder cancer. Using the multi-scan imaging and image segmentation 
techniques described above, the region of the bladder wall can be readily delineated. 
Regions of normal bladder tissue generally exhibit a substantially uniform texture 
feature. However, if a tumor is present in the region, the uniform texture feature will 

5 be interrupted. Thus, using texture analysis to evaluate the wall of the bladder, a 

region which may exhibit a tumor will present itself as a disturbance, or "noisy region" 
within the uniform texture. 

The texture of a region can be represented by a probability distribution 
function (PDF) which characterizes the intensity correlation between voxels within a 

0 defined range. A two-dimensional PDF can be used to represent a texture feature. 
Such a PDF characterizes the correlation between two closest voxels along all 
directions. To estimate the PDF, the intensities of any two closest neighbor voxels in a 
region of interest can be recorded as a sample vector for the region of interest (e.g., 
context). Using a number of such sample vectors, a cumulating distribution function 

5 (CDF) can be generated which estimates the PDF for that context. For each voxel, 
sample vectors within a range of its neighbor can also be used to generate a local CDF. 

A statistical test, such as a Kolmogorov-Smimov test, can be applied to 
the CDF to determine whether the CDF of the context and the local CDF are 
statistically equivalent, e.g., within a predefined confidence level. If so, the local 

0 texture feature around the current voxel is regarded as identical to the context. 
Otherwise, the current voxel exhibits a different texture feature firom that of the 
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context and may be regarded as a potential abnormality, such as a tumor. The level of 
confidence used to detemiine whether a voxel is statistically equivalent to the context 
can be varied to increase or decrease the sensitivity of detection. 

In an altemative method of applying the PDF or CDF for texture 
5 analyis, each CDF or PDF can be regarded as a point in a functional linear space. The 
distance between two CDF's or PDF's in that space can be measured, such as in terms 
of the Skorohold metric. This distance provides*a measure of the degree of similarity 
of PDF's. For example, the distance between a local CDF and the context CDF can be 
calculated and the resulting distance can be compared to one or more distancing 

1 0 thresholds. If the distance is large, the local texture may be considered different from 
the context, which can indicate that such a voxel belongs to a region with a potential 
abnonnality or tumor. Preferably, the distancing thresholds are determined based on 
evaluation of a statistically sufScient known data sets. 

The distance calculated above can be used with visualization techniques 

15 and volume rendering techniques, such as those described herein. For example, a 

feature volimie dataset having a size comparable to the original dataset can be created. 
The intensity for each voxel in the new dataset can then be assigned based upon the 
distance between the local CDF and the CDF of the context. When this three 
dimensional volume dataset is viewed . through volume rendering techniques, the 

20 regions which contain suspected tumors will exhibit a higher image mtensity than the 
surrounding area. 

As was discussed above in connection with automatic detection of 
polyps, the surface of a lumen can be represented as a continuously second 
differentiable surface in three dimensional Euclidean space, such as by using a C-2 

25 smoothness surface model. In such a model, each voxel on the surface of the colon has 
an associated geometrical feature which has a Gauss curvature, referred to as Gauss 
curvature fields. For various organs, certain expected local features can be 
characterized by distinct curvature templates. For example, in the context of the colon, 
the ejected local features include smooth curve surfaces, ring folds, convex hills from 

30 a smooth surface and plateaus from a smooth surface. These last two local features 
may be indicative of a polyp or tumor. Accordingly, by searching the Gauss curvature 
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fields for specific predetennined local feature templates, polyps, tumors and other 
abnormalities of interest can be automatically detected. This use of surface geometry 
to perform computer assisted diagnosis can be used alone, or in conjunction with the 
texture-based CAD techniques described above, 
5 As an alternative to rendering and viewing the entire organ or region of 

interest, the surface under observation can be partitioned into small areas, or patches, 
which are defined by the local curvature templates. Each patch should contain voxels 
which have a common geometry feature, or curvature template. A single viewing point 
is then determined for the patch which allows all voxels of the patch to be observed. 

10 The patches are then assigned a priority score indicating the probability that the patch 
' represents a polyp or other abnormality. The patches can then be observed 
individually, in priority order, rather than requiring the operator to navigate the entire 
organ volimie to search out suspect areas. Of course, a preferred diagnostic system 
includes the ability to toggle between views such that an operator can readily change 

1 5 from viewing a patch to viewing the patch in the context of the organ. Alternatively, 
these two views can be presented simultaneously. Again, the texture based approaches 
can be used to supplement this approach. By mapping the results of texture analysis 
onto the patches being observed, the texture information can also be observed and used 
in diagnoses. 

20 The foregoing merely illustrates the principles of the present imaging 

and examination systems and methods. It will thus be appreciated that those skilled in 
the art will be able to devise numerous systems, apparatus and methods which, 
although not explicitiy shown or described herein, embody the principles of the 
invention and are thus within the spirit and scope of the invention as defined by its 

25 claims. 

For example, the methods and systems described herein could be 
applied to virtually examine an animal, fish or inanimate object. Besides the stated 
uses in the medical field, applications of the technique could be used to detect the 
contents of sealed objects which caimot be opened. The techniques can also be used 
30 inside an architectural structure such as a building or cavem and enable the operator to 
navigate through the structure. 
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CLAIMS 

1 . A method for performing virtual examination of an object comprising: 

performing at least one scan of an object with the object distended by 
the presence of a contrast agent; 
5 performing at least one scan of the object with the object relieved of the 

contrast agent; 

converting the scans to corresponding volume datasets comprising a 
plurality of voxels; 

performing image segmentation to classify the voxels of each scan into 
10 a plurality of categories; 

registering the voliune datasets of each scan to a common coordinate 

system; 

displaying at least two of the volume datasets in a substantially 
simultaneous manner; and 
1 5 performing virtual navigation operations in one of the volume datasets 

and having the corresponding navigation operations take place in at least one other 
volume dataset 

2. The method for performing virtual examination according to claim 1 , wherein 
the at least one scan of the distended object includes a transverse scan and a coronal 

* 

20 scan of the object 

3. The method for performing virtual examination according to claim 2, wherein 
the at least one scan of the relieved object includes a transverse scan and a coronal scan 
of the object. 

4. The method for performing virtual examination according to claim 3, wherein 
25 the object is a bladder. 
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5. The method of performing virtual examination according to claim 4, wherein 
the scans are computed tomography scans. 

6. The method of performing virtual examination according to claim 4, wherein 
the scans are ultrasound imaging scans. 

5 7. The method of performing virtual examination according to claim 4, wherein 
the scans are magnetic resonance imaging scans. 

8 The method of performing virtual examination according to claim 7, wherein 
the contrast agent is urine. 

9. The method for performing virtual examination according to claim 1, wherein 
10 the at least one scan of the relieved object includes a transverse scan and a coronal scan 

of the object 

1 0. The method for performing virtual examination according to claim 1 , wherein 
the object is a bladder. 

1 1 . The method of performing virtual examination according to claim 1 0, wherein 
15 the scans are computed tomography scans. 

12. The method of performing virtual examination according to claim 1 0, wherein 
the scans are ultrasound imaging scans. 

13. The method of performing virtual examination according to claim 10, wherein 
the scans are magnetic resonance imaging scans 

20 14. The method of performing virtual examination according to claim 13, wherein 
the contrast agent is urine. 
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15. The method of perfonning virtual examination according to claim 1, further 
comprising evaluating the at least one scan with the object distended and the at least 
one scan with the object relieved to identify regions where contrast is more visible in 
one of said scans and evaluating the scan with more contrast in a region of interest to 

5 determine physiological characteristics of the object 

16. The method of performing virtual examination according to claim 1 5, wherein 
said step of image segmentation includes classifying voxels based on local intensity 
vectors of the voxels. 

1 7. The method of performing virtual examination according to claim 1 6, wherein 
10 the step of image segmentation further includes usmg a region growing algorithm to 

identify regions of the object based on the classified voxels. 

1 8. The method of performing virtual examination according to claim 1, further 
comprising partitioning the volume image datasets into a plurality of regions related to 

« 

the coordinate system. 

15 19. The method of performing virtual examination according to claim 1 8, wherein 
the plurality of regions include eight regions defined in a three dimensional coordinate 
system. 

20. A method for perfonning virtual examination of an object comprising: 
20 performing an imaging scan of the object to acquire image scan data; 

converting the acquired image scan data to a plurality of voxels; 
interpolating between the voxels to generate an e}q)anded dataset; 
performing image segmentation to classify the voxels into a plurality of 

categories; 

25 extracting a volume of the object interior from the expanded dataset; 

generating a reduced resolution dataset from the expanded dataset; 
storing the expanded dataset in a tree data structure; 
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rendering images for the expanded dataset and reduced resolution 

dataset; and 

selecting at least one of the reduced resolution dataset or expanded 
dataset renderings for display. 



5 21 . The method for performing virtual examination of an object of claim 20, 

wherein the selecting step comprises: 

selecting the reduced resolution dataset during image interaction; and 
selecting the expanded dataset rendering if no image interaction has 

occxirred in a predetermined time period. 



10 22. The method for performmg virtual examination of an object of claim 20, 
wherem the imaging scan is a computed tomography scan. 

23 . The method for performing virtual examination of an object of claim 20, 
wherein the imaging scan is a magnetic resonance imaging scan. 

24. The method for performing virtual examination of an object of claim 20, 
15 wherein the imaging scan is an ultrasoxmd imaging scan. 



25. The method for performing virtual examination of an object of claim 20, 
wherein the object is the larynx. 

26. The method for performing vutual examination of an object of claim 20, 
wherein the tree structure is a binary space partition tree structure. 

20 27. A method of performing virtual angiography comprising: 

acquiring imaging scan data including at least a portion of the aorta; 
converting the imaging scan data to a volume representation includmg 
plurality of voxels; 
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segmenting the volume representation to classify the voxels into one of 
a plurality of categories; 

analyzing the segmented volume representation to identify voxels 
indicative of at least a portion of an aneurysm in the aortic wall; and 
S generating at least one closing surface around the voxels indicative of at 

least a portion of an aneurysm to estimate the contour of the aneurysm, 

28. The method of perfomiing virtual angiography of claim 27, wherein the 
imaging scan is a computed tomography scan. 

29. The method of perfomiing virtual angiography of claim 27, wherein the 
1 0 imaging scan is a magnetic resonance imaging scan. 

30. The method of performing virtual angiography of claun 27, wherein the 
segmenting operation classifies voxels in at least the categories of blood, tissue, and 
calcium deposits. 

3 1 . The method of performing virtual angiography of claim 27, further comprising 
1 5 estimating the volume of the aneurysm using the generated closing surfaces. 

32. The method of performing virtual angiography of claim 27, further comprising 
generating a navigation path through the aortic lumen. 

■ 

33 . The method of performing vntual angiography of claim 32, further comprismg 
20 estimating the lei^th of the aneurysm based on the navigation path. 

34. A method of performing virtual endoscopy of a blood vessel comprising: 

acquiring imaging scan data including at least a portion of the vessel; 
converting the imaging scan data to a volume representation including a 
plurality of voxels; 
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segmenting the volume representation to classify the voxels into one of 
a plurality of categories including the categories of bloody tissue, and calcium deposits; 
and 

generating a navigation path through the vessel 

5 35. The method of performing virtual endoscopy of claim 34, wherein the vessel is 
a carotid artery. 

36. The method of performing virtual endoscopy of claim 34, further comprising 
the step of determining the diameter of the carotid artery along the navigation path to 
identify regions of narrowing. 

10 37. The method of performing virtual angiography of claim 34, wherein the 
imaging scan is a computed tomography scan. 

38. The method of performing virtual angiography of claim 34, wherein the 
imaging scan is a magnetic resonance imaging scan. 

39. A method of determining the characteristics of a stent graft using virtual 
15 angioscopy, comprising: 

acquiring imaging scan data including at least a portion of the aorta; 

« 

converting the imaging scan data to a volume representation including a 
plurality of voxels; 

segmenting the volume representation to classify the voxels into one of 
20 a plurality of categories; 

analyzing the segmented volume representation to identify voxels 
indicative of at least a portion of an aneurysm in the aortic wall; 

generating at least one closing surface around the voxels indicative of at 
least a portion of an aneurysm to estimate the contour of the aneurysm; 
25 identifying the location of the endpoints of the aneurysm contour; 
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calculating the length between the endpoints of the aneurysm contour to 
determine the length of the stent graft; and 

calculating the diameter of the aortic lumen at the endpoints of the 
aneurysm contour to determine the required outside diameters of the stent graft. 

5 40. The method of determining the characteristics of a stent graft of claim 39, 
further comprising determinii^ the angle of interface of the aneurysm and normal 
aortic lumen to determine an angular direction of a corresponding end of the stent 
graft. 

41 . The method of determining the characteristics of a stent graft of claim 39, 
10 further comprising locating arterial branches proximate the aneurysm to determine a 

maximum length of the stent graft. 

42. The method of determining the characteristics of a stent graft of claim 41, 
wherein the arterial branches proximate the aneurysm include at least one of the renal 
and femoral arterial branches. 

15 43. The method of determining the characteristics of a stent graft of claim 39, 

further comprising conducting a virtual biopsy of the aortic region proximate the ends 
of the aneurysm to determine the nature of the tissue at the anticipated graft interface 
locations. 

44. A method of defining a skeleton for a three dimensional image representation 
20 of a hollow object formed with a plurality of voxels comprising: 

identifying a root voxel within the hollow object; 

generating a distance map for all voxels within the hollow object, the distance 
map being formed using a 26-connected cubic plate of neighboring voxels having 
Euclidian weighted distances; 
25 identifying voxels having a local maxima in the distance map as endpoints of 

branches in the hollow object; and 
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for each local maxima voxel, determining a shortest connected path to one of 
the root voxel or a previously defined shortest path. 

45 . The method of defining a skeleton for a tbree dimensional image representation 
of claim 44 further comprising performing multi-resolution data reduction to the three 

S dimensional image representation to generate a reduced data set for the generating and 
identifying operations. 

46. The method of defining a skeleton for a three dimensional image 
representation of claim 44 fiirther comprising centralizing the shortest paths within the 
respective branches of the object. 

1 0 47. The method of defining a skeleton for a three dimensional image 

representation of claim 44, wherein the object includes at least one blood vessel. 

48. The method of defining a skeleton for a three dimensional image representation 
of claim 44, wherein the object includes the airways of a lung. 

49. The method of defining a skeleton for a three dimensional image representation 
15 of claim 44, wherein the object includes the bladder. 

50. The method of defining a skeleton for a three dimensional image representation 
of claim 44, wherein the object includes the spinal cord of a vertebrate animal. 

51. A method of performing computed assisted diagnosis of a region of interest, 
comprising: 

20 acquiring imaging scan data including at least a portion of the region of 

raterest; 
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converting the imaging scan data to a volume representation including a 
plurality of voxels, at least a portion of the voxels representing a surface of the region 
of interest; and 

analyzing said portion of voxels representing a surface for at least one 
5 of a geometric feature and a textural feature indicative of an abnormality. 



52. The method of performing computed assisted diagnosis according to claim 5 1 , 
wherein the textural feature is included in a probability density function characterizing 
a correlation between two voxels of the portion of voxels. 



53 The method of performing computed assisted diagnosis according to claim 52, 
10 wherein the two voxels are adjacent voxels. 



54. The method of performing computer assisted diagnosis according to claim 52, 
wherein intensities of said portion of voxels are used to generate an estimate of the 
probability density function. 



55. The method of performing computer assisted diagnosis according to claim 54, 
1 5 wherein a plurality of voxel intensities are used to generate a cumulating distribution 

function of the region of interest and a local cumulating distribution Amction, and 
wherein the local cumulating distribution Auiction is compared against the context 
cumulating distribution function to identify regions of abnormality. 

56. The method of performing computer assisted diagnosis according to claim 55, 
20 wherein a distance is determined between said local cumulating distribution function 

and said context cumulating distribution function, the distance providing a measure of 
abnormality. 



57. The method of performing computer assisted diagnosis according to claim 56, . 
wherein the distance is used to assign intensity values to the voxels representing a 
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surface of the region of interest and wherein said method further comprises displaying 
said voxels such that variations in intensity rq)resent regions of abnormality. 

58. The method of performing computer assisted diagnosis according to claim 57, 
wherein the region of interest includes the colon and wherein the abnormality includes 

5 polyps. 

59. The method of performiiig computer assisted diagnosis according to claim 51, 
wherein the region of interest includes the aorta and wherein the abnormality includes 

» 

abdominal aortic aneurysms. 

60. The method of performing computer assisted diagnosis according to claim 51 
10 wherein the surface is represented as a second diflferentiable surface where each surface 

volume imit has an associate Gauss curvature and wherein said Gauss curvatures 
combine to form said geometric features. 

■ 

6 1 . The method of performing computer assisted diagnosis according to claim 59 
wherein a plurality of predetermined geometrical feature templates are defined and 

15 wherein the geometric features of said surface are compared to said templates to 
determine a geometric feature classification. 
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